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Two   problems   are  presented  which  are   linear  on  two   adjacent    inter- 
vals  but   not  on  their  union.      These  problems   are  associated  with   the 

(AX  +  BF,    0  <t  <ti 
differential   equation  X  =    7qX  +  DF     t  <t<T    ■   wnere  x   *-s   fc^e  matrix 

(    J,   where   F   is   a   2   x   1  matrix,    and   where  A,B,C,    and   D  are   2x2  matrices 

of   functions  of  t.      t,    is   a  variable,   hence   the  differential    equation   is 

non-linear.      Problems   associated  with  this  differential    equation   are 

called   semi-linear. 

In  the  first   problem,    a   condition  is   found  on  t,    and   F  which  must 
be  satisfied  whenever  x(T)    is   to   be  a  maximum  with  y(T)    fixed.      In  the 
second   problem,    conditions  on  F  and   t^   are   found   which  must   be  satisfied 
for  x(T)    to   be  a  maximum  for   a   fixed  y(T)    and    for  a   fixed   x(t^) .      A  numer- 
ical  routine   is   developed  which  yields    successive  approximations   to   the 
maximum  value  of  x(T). 

The  basic   theory  of  the  methods    is   presented,    and   the   problems   are 
developed   in  the   context  of  optimum  control. 
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I.      Introduction. 

Nonlinearities  make  the   solution  of  optimum  control   problems  more 
difficult.      In  general,    different   techniques  must  be  developed   for  each 
new  nonlinear  problem. 

The  particular  problems  to   be  discussed   are  those   in  which  a  set  of 

functions   x.   of  a  variable  t   is   related  to   a   second   set   of   functions   f . 
1  1 

of  t   through  a  differential    equation  which   is   linear  on  each  of  the    in- 
tervals   (0,t-i)    and    (t-pT).      If  t-.    is   a  variable,   the  differential    equa- 
tion,  though   linear  on  each  of  the   intervals   (0,t-,)    and    (t-,,T),    is 
non-linear  on  the  interval    (0,T).      Such  a  differential    equation  will   be 
called   a  semi-linear  equation;    problems   in  which  a   semi-linear  differen- 
tial  equation  occurs   will   be   called   semi-linear  prdblems.      Such  problems 
arise   in  the  theory  of  optimum  control;    the   problems  to  be  discussed 
will  be   in  the  terminology  of  optimum  control   theory. 

The  specific   semi-linear  problems  to  be  considered   arise  from  the 

.  fAX    +    RF      O'C  tc  t 

differential    equation  X   =   ) „  r  *  *    ,   where  X,A,B,C,D,    and   F 

M  (CX  +   DF,    t1<t<T 

are  described  as   follows:      X  is   the  matrix    (~1   whose  elements  are  con- 
tinuous   functions  of   t;   X   is  the  matrix   I* J  of  derivatives  of  x  and  y 
with  respect  to   t.      A,B,C,    and  D   are  2x2  matrices  whose   elements   are 
functions  of  t  which   are  piecewise  continuous   and   bounded  on  the   inter- 
val   (0,T).      In  addition   it   is   necessary  that  B   and  D   be  nonsingular.      F 

is   the  matrix  lf  J    of   functions  of  t;    associated  with   F  there   is   some 
2 

constraint,  specified  in  each  problem.  F  is  called  the  matrix  of  con- 
trol variables,  and  those  matrices  F  meeting  the  specified  constraints 
are  called   allowable. 

Associated  with  the  above  differential   equation  are  certain  boundary 
conditions   which  X(t)    is  to   satisfy.      Generally  the  initial   point    is 


given.   In  addition,  one  or  more  components  of  X  may  be  specified  at  some 
value  of  t  on  the  interval  (0,T).   Curves  satisfying  the  given  boundary 
conditions  and  for  which  F  is  allowable,  are  called  admissible. 

In  the  terms  that  have  been  defined  above,  and  for  the  above  differ- 
ential equation,  the  semi-linear  problems  to  be  solved  are  the  following: 
(1)  Find  conditions  on  F  and  on  t-,    which  must  be  satisfied  if  x(T)  is  to 
be  a  maximum  for  a  given  y(T),  where  t,  must  be  determined  and  where  x(t,) 
and  y(t,)  are  unspecified.   (2)  Add  the  constraint  tV^-t  x(t  )  is  to  be 
specified,  for  a  value  of  t.  to  be  determined.   Then  (a)  find  conditions 
on  F  and  t-,  which  must  be  satisfied  for  a  curve  to  be  admissible;  (b)  find 
additional  conditions  on  F  and  on  t,  such  that  an  admissible  curve  will 
yield  the  maximum  x(T) ;  and,  (c)  develop  a  method  of  successive  approxima- 
tions which  will  give  values  of  x(T)  converging  to  this  maximum. 

The  following  results  are  obtained:   In  the  first  problem,  a  neces- 
sary condition  at  t-,  is  derived.   In  the  second  problem,  conditions  for 
a  maximum  are  derived,  and  a  numerical  method  is  given  for  obtaining  this 
maximum.   These  results  are  another  step  in  the  solution  of  nonlinear  prob- 
lems of  optimum  control. 

I  wish  to  thank  Professor  Fpulkner  for  the  encouragement  and  help  he 
has  given  me  and  for  the  permission  to  use  class  notes  of  his  course, 
Differential  Equations  of  Optimum  Control. 


II.   General  terminology  and  some  sufficient  conditions  for  the  linear 
problem. 

In  this  section  we  define  the  linear  problem  and  explain  the  termi- 
nology to  be  used.  We  then  derive  some  sufficient  conditions  for  the 
linear  problem. 

Let  us  consider  the  differential  equation 

(1)  X~  AX  +  BF 

where  X  is  the  n  x  1  matrix  (x.)>  where  A  and  B  are  the  n  x  n  matrices 
(a-  •)  and  (b.  .)  ,  respectively,  and  where  F  is  the  n  x  1  matrix  (f.).   We 
will  always  assume  that  B  is  nonsingular,  although  the  results  apply  to 
many  problems  where  B  is  singular.   The  a. .,  the  b. .,  and  the  f-    are 
functions  of  t  that  are  bounded  and  piecewise  continuous  on  the  inter- 
val (0,  T) .   The  a-  -  and  the  b.  .  are  given  .functions;  the  f.  are  func- 
tions satisfying  a  given  constraint  but  otherwise  unspecified.   The  f. 
are  often  called  control  variables  and  the  x-  state,  or  dependent 
variables. 

Associated  with  the  differential  equation  (1)  are  certain  boundary 
conditions  which  we  want  X(t)  to  satisfy.   Generally  the  initial  point 
is  given.   In  addition  we  may  specify  that  certain  elements  of  X  take  on 
stated  values  at  various  fixed  points  of  the  interval  (0,  T) .   Note  that 
the  differential  equation  is  linear,  since  A  and  B  are  functions  of  t 
only. 

If  we  multiply  both  sides  of  equation  (1)  on  the  left  by  the  1  x  n 

T  i 

matrix  K  (the  transpose  of  the  matrix  K) ,  whose  elements  k  are  func- 
tions which  have  not  yet  been  specified,  and  integrate  from  t=0  to  t=T, 

we  get 

T 

(2)  \  KAX  dt  =  NK'-AX  dt  +  \KTBF  dt 


\  KTX  dt  =  \KTAX  dt  +  W1 


Integrating  the  left  side  by  parts  and  rearranging  terms,  we  get 
T    T  T 


( 3)  KTX 


=     \  (K      +    KTA)X  dt    +     J\KTBF  dt. 


Now  let  us  choose  K     =   -   K  A  and  take  the  transpose;    we  get 

(4)  K  =   -ATK. 

This    is  the  adjoint   equation  associated  with  equation   (1).      Choosing   K 
as   a   solution  to   the  adjoint   equation  and   substituting   into   equation    (3), 
we  get 


(5)  KTX 


T        T 

=   Q   KTBF   dt. 


0        0 


1  IT 

We  can  choose  a   solution  K     to   the  adjoint   equation  so   that    K     (T)    is 

the  matrix   (1     0     0      ...    0).      If   we  substitute  this   solution   into   equa- 
tion   (5)    and   rearrange  terms,    we   get   a  solution  for  x-,(T),    namely 

T 
(6)      xx(T)  =  K1T(0)X(0)  +  J  KlT(t)BF  dt. 

0 

1  9     "3  n 

In  the  same  way  we  chose  K  ,  we  can  choose  K  ,  K  ,  ...  ,  K  such  that 

iT 
K   (T)  =  (b^S^.  •••  6  • )  ,  where  6.  .  equals  zero,  if  i/  j ,  or  one,  if 

i=j.   Suppose  now  that  we  have  the  n  linearly  independent  column  mat- 
rices K1  chosen  above.   Such  a  set  of  n  linearly  independent  matrices 
is  called  a  fundamental  set  of  solutions  for  equation  (4).   If  we  form 
the  n  x  n  matrix  /\   whose  i'th  column  is  the  matrix  K  ,  then  we  may 
write  every  solution  to  equation  (4)  in  the  form  /CC,  where  C  is  an 
n  x  1  matrix  of  constants.   Conversely,  every  product  of  the  form/^C, 
being  a  linear  combination  of  solutions  to  equation  (4),  is  also  a 
solution.   Furthermore  fK.   is  itself  a  solution.   This  is  shown  as  fol- 
lows:  Take  any  product  y^C,  where  C  is  an  n  x  1  matrix  of  constants 
and  where  /\,   is  the  matrix  defined  above.   This  product  is  a  solution 


to  equation  (4).   Hence /\,C  =  -A  f^C.      But  this  equation  is  valid  for 
every  C  and  hence  for  the  C  whose  transpose  is  (1  0  0  ...0).   Substitu- 
ting this  C  into  the  above  equation,  we  see  that  the  first  column  of  rC, 
is  the  same  as  the  first  column  of  -A  /{.      But  we  could  have  equally  well 
chosen  the  i'th  element  of  C  as  one  with  the  others  zero;  we  hence  could 
have  found  that  the  i'th  columns  of  both  sides  were  equal,  for  isl,  2,  ..., 

s 

n.      Hence /\=  -A  /^,    i.e.  /t   is   a   solution  to   the  adjoint   equation,   which 
is  what  we  wanted  to   show.      Furthermore,    if    6  JC  is   a  small   variation   in 
*]/  ,    and   if  6  >T  is   the  corresponding  variation   in  /^,    then   J£  +  b  /\  = 
-A  (/L  +  6ft),    i.e.     6/^=  -A   6^,    since    /^  =   -A    f\.      Hence  the  varia- 
tions   in  J\,  are  related  to   the  variations    in  /^  by  the  adjoint  equation. 
We  have  shown  that   solutions   to   the   adjoint  equation  can  be  used   to 
find  solutions  to   the  equation  X  =  AX  +  BF  when  F  is  known.      The  next 
problem  we  are  concerned  with   is   that  of   finding   F  so   that   for  a  given 
X(0)    and   a  given  T,    a  specified   component  of  X(T)    is   a  maximum.      F  can 
be  regarded  as   a  column  matrix  of   forcing  functions;    these  forcing  func- 
tions  are  the  components  of   the   forcing  function  vector.      There  may  be 
one  of   several  types  of   constraints   on  F.      In  rocket   thrust   scheduling, 
for   example,    the   acceleration  possible  at   any  one  time   is   limited  by  a 
function  of   the  mass  of   fuel  on  board  at   that   time.      If  we  call  this 
function    <p(t) ,    the  corresponding   constraint  on  F   is   that    |Fp<P(t)    , 
where   Jf|    is  the  square  root  of  the  sum  of  the  squares  of  the  elements 
of  F.      Another  type  of   constraint  on  F   is    IfJ^la    I,    where  a.    is   a 
given  function  of  t.      Problems  of  this   second  kind   are  called  bang- 
bang  control   problems ,  £ij;i      In  every  problem   it   is   necessary  to   state 
the  constraint  on  F;    forcing  functions   which   satisfy   the  stated   con- 
straints will   be  called   allowable.      Solution  curves   to   equation    (1) 
for  which  F  is   allowable  will   be  called   allowable  curves:    allowable 


curves  satisfying  the  differential  equation  and  satisfying  the  given 
boundary  conditions  will  be  called  admissible. 

We  want  a  method  for  choosing  an  allowable  F  such  that  X(0)  is  a 
given  point  and  such  that  a  specified  element  of  X(T)  is  a  maximum.   A 
most  important  principle  enables  us  to  state  conditions  which  F  must 
satisfy  whenever  the  desired  maximization  takes  place,  namely  Pontry- 
agin's  Maximum  Principle;   The  control  variables  must  be  chosen  from  the 
set  of  allowable  controls  so  as  to  maximize  a  scalar  product  of  some 
solution  to  the  adjoint  equation  and  the  forcing  function  vector  at  every 
time  t.-[2j$  We  will  use  this  principle  extensively  in  the  following 
pages. 

Consider  the  following  problem  and  see  how  the  maximum  principle 

may  apply  to  it.   Suppose  that  we  want  a  curve  beginning  at  some  fixed 

point  at  time  t  =  0  on  which  some  element  of  X  at  t=T,  say  x,(T),  is  a 

maximum  and  such  that  on  the  curve  these  three  conditions  hold:   First, 

F  is  allowable.   Second,  the  curve  satisfies  the  differential  equation 

X  =  AX  +  BF  for  all  values  of  t  between  t=0  and  t=T  where  T>  0  is  given. 

Third,  x^(T)  =  x-T,  for  i=2 ,  ...  ,  N,  where  x.   are  given,  and  where 

KN^n,  and  where  x.,^,  (T) ,  ...  ,  x  (T)  are  unspecified.   N=l  means  that 

jn+1  n 

no  values  of  the  x^  are  given  at  t=T;  N=n  means  that  all  values  except 
X-.  are  given  at  t=T.   A  proof  given  by  Faulkner  ^Sj,  proves  that  the  fol- 
lowing hypotheses  are  sufficient  for  a  curve  C*  to  yield  a  maximum  x,(T). 
Suppose  that  we  have  found  a  curve  C*  with  forcing  function  F  and  have 
found  at  the  same  time  a  solution  K*  to  the  adjoint  equation  which  to- 
gether satisfy  the  following  hypotheses: 

HI.   C   is  admissible:  it  begins  at  XQ  and  ends  on  the  manifold  de- 
fined by  X2(T)  =  x2„,  ...  ,  xN(T)  =  x  _,,  and  F*  is  allowable. 
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H2.      F*  maximizes   K*TBF  for  all   allowable   F,    i.e.    K*TBF*is  K*TBF  for 
all  allowable  F. 

H3.      k^(T)>0,   and  k£(T)    =  0   for   i>  N.      No   restriction   is   put  on 
k?(T)   for   i  =   2,    ...    ,N. 

■k 

THEOREM;   C     furnishes  the  desired  maximum  of  x^(T) . 

T   T 

T 

Proof:  We  have  shown  that  K'-X 


=  \K  BF  dt  for  every  solution  K  to  the 
0   0 


•k 

adjoint   equation.     Hence  for  the  particular  solution  K     and  the  matrix 
X*  we  have,   on  the  path  C    , 

(7)  ki(T)x*(T)    +  k*(T)x2T   +    ...    +  k*(T)xNT  =    (K*TX*)T 

T 

*t  *       r  *t 

=    (K     X  ).    +  \K     BF  dt. 
°        0 

Consider  any  other  admissible  path  C*  with  F  =  F  .  For  this  path  and 


for  K  =  K*  we  have 


(8)  k^(T)x|(T)  +  ...  +  kN(T)xNT  =  (K  X  )Q  +  ^K  BF  dt. 

0 

Subtracting  equation  (8)    from  equation   (7),   we  get 

T 

(9)  k^(x1-x1)   T|\(K     BF   -K     BF   )      dt  £    0. 

0 

Hence,  since  k1(T)>  0,  x-^CT)^  xj  at  t  =  T.   Hence  the  given  hypotheses 

JL 

are  indeed   sufficient   to  give  a  maximum  x,(T).      Note  that   since  k-(T)    = 

•k 

0,    for   j=N+l,    ...    ,    n,    the  unspecified  elements  of  X   (T)    play  no  part   in 
the  solution. 

Having  developed   a  useful   sufficiency  condition,    let   us   now  consid- 
er the  case  where  X   is  the  matrix    (  J,   where  A  and  B  are   2x2  matrices, 

where  F    is   the  matrix  If),    and  where  we  want   to  maximize  x(T)    subject   to 

2  2 

the  following  constraints:      First,    (f^)      +    (f2)      =   1.      Second,   T>0   is 


11 


fixed.      Third,    solution  curves  must   start   at   X     and   end  on  the   line  y(T)    : 

Vf,   where  yf   is  given.      For   this   problem  the  admissible  curves   are  allow- 

« 
able  curves   satisfying  the   third   constraint. 

1  2 

Let   us   choose  the  solutions   K      and  K     to   the   adjoint   equation  having 

K  (T)    =  (q)    and   K^(T)    =    (J.      Using   these  solutions   in  equation   (6),    we 

T 
IT  f    IT 

(10)  x(T)    =    (K     X)Q    +  \K     BF  dt 

0 


and 

T 

2T  f    2T 

(11)  y(T)    =    (K     X)Q   +  ^K     BF   dt. 

0 

2  2 

Since  the   given  constraint  on  F  is   that    (f^)    +(f2)      =   1 ,    we  may  choose 

fi    =  cos   9,    f~   =   sin  9,   where  9-^is    a   function  of   t.      With   this   substitu- 
tion,   a  variation   in  0  will   cause   a  variation   in   F  which  will    in  turn 
cause  a  variation   in  x  and    in  y.      Let   us   call   the  variation   in  0,  60  and 
the  variations   in  x  and   in  y   at   t=T,  fix  and  6y,    respectively.      We  then 
get,    from  equations    (10)    and    (11), 

T 


(12)  6x  =  Jk     Si  69  dt 


and 


f    2T     i  F 

(13)  6y  =  JK     B-J|-6e  dt, 

0 

for  sufficiently  small  values  of   &0.      Let   us   now  choose  a   special  varia- 

1T    3  F  9T    d  F 

tion  of  the   form     0  =   e^K     B  X g      +   e^K^B  S  g     where  e,    and  e     are  con- 
stants yet   to  be  chosen.      Substituting  this  variation   in   equation    (12), 
we  get  T 

(14)  6x  =   \(K1TB-||-)(e1K1Trj|-  +  e2K2TB-|-|-)dt 

0 


12 


and  a  similar  expression  for  &y.   These  can  be  simplified  if  we  let 

T 


iij=Kr4f-><KJT4f>dt 


o 

Using  this   substitution  we  get 

(15)  6x  =  e1I11    +  e2I12 

and 

21  99 

(16)  6y  =  exI        +  e2I      . 

12  21  22 

Note  first  that    I        =   I  Note  next   the  consequences    if   I      =0.      If 

92  2T    <&  F  91  c 

I        =0,   then  K     B  must  be  identically   zero;    hence  I        and  by  are 

both   identically   zero;    hence  y(T)    is   not   affected  by  a  change   in  the  con- 
trol variable  0.      We  describe  such   a  situation  by   saying  that   the  curve 
furnishes    a   stationary  value  for  y(T).      We  will   not   want  y  to   have  a 

stationary  value,    and  we  shall  henceforth  assume  for  our  curve  that   I    ¥ 

9T     ^  F 
0,   or,    equivalent  ly,   that   KB  is   not    identically   zero  on  the   inter- 

val    (0,T). 

Since  T  is   fixed,    equations   (15)    and    (16)    give  the  total   differen- 
tials of  x  and  of  y.      If  we  replace  6 y  by  y^  -  y(T) ,   we  can  generally 
find  values   for  e-,    and  e„   which  give  a  variation  which  will  make  the  re- 
sulting curve  admissible. 

Suppose  now  that  we  have  an  admissible  curve  and  that  we  want   to 
find  conditions  on  the   I        such  that   x(T)    is   a  maximum.      For  an  admissi- 
ble  curve  we  have 

(17)  &x  =  e-jj        +  e2I 


and 


(18)  0   =   exI21    +  e2I22. 


Equations    (17)    and    (18)   can  be   solved   for  £x>0  only   if  the  rank  of 


13 


'6x  I11  I12' 

1  is  equal  to  the  rank  of  /        I.   But  the  rank  of  the  first 


0  I21  I22 


/bx     0   0  \ 
(  0   I21  I22] 


matrix  is  the  same  as  the  rank  of  |      01   o0  |  .   The  rank  of  this  last 

21   22 
matrix  is  one  greater  than  the  rank  of  the  matrix  (I   I  );  hence  equa- 
tions (17)  and  (18)  cannot  be  solved  for  <5x  whenever  the  rank  of 

IU  IU\  21   22 

I  is  equal  to  the  rank  of  (I   I  ).   But  this  will  be  true  only 

j-21   I22/ 

when  the  determinant    II      I    is    zero.      Hence   if  the  admissible  curve  yields 
a  maximum   for  x(T),    then      I    J      =  0. 

Let   us   consider  some  of  the  consequences  of   this   condition.      If  the 
determinant   |lL^|    ■  0,    then  the  matrix   (I1^)    has   a   zero  eigenvalue.      Let 

(c  i  \  £  -j 

)  is   an   eigenvector  of    (I   J) 

corresponding  to   this    zero   eigenvalue.      Then  let   us   consider 

d9)    j  =  W1tb4£-  ♦  c X"*4&2  at. 

Clearly  J£0.   But  J  =  (c1)2I11  +  2c1c2I12  +  (c2)2I22 


■ ""''  ('» $  fc) 

C) 


ButJ  \   is  an  eigenvector  corresponding  to  the  zero  eigenvalue  of  (I  ), 

hence  J  =  0.   Therefore  c,K  B-r-s—  +  c„K  B  &  F  =  0.   This  equation 

1^8     2^0 

appears  frequently  in  the  literature  on  calculus  of  variations  -  it  is 

known  as  the  Euler  equation.   Curves  whereon  the  Euler  equation  is  sat- 

i         *      1      2 
isfied  are  frequently  called  extremals.   Hence  K  =  c^K  +  c2K  is  a 

This  is  the  definition  given  by  Bolza  j**]  .   Another  definition  is 
given  by  Bliss  [5] 
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solution  to   the  adjoint   equation  which  gives   an  extremal. 

We  saw  that   if  c,  >  0   and   if  9  =  0     maximizes   K     BF  then  xft1)    is   a 
maximum.      The  condition  K     B  g  g     =  0    is   a  necessary  condition  for  x(T) 
to  be  a  maximum.     We  can  and  will   choose  c-,>  0.     We  will   also   see  that 
it   is   necessary  that   8  maximize   K     BF;    this   is   the  Weierstrass   condition. 
In  this  problem  the  Weierstrass   condition  is  necessary  and   sufficient 
for  an  admissible  curve  to   furnish  the  desired  maximum. 

Let   us   assume  that  the  rank  of   (I   J)    is   at   least  one  for  all  curves 
which  arise:   this   is  called  normality  for  the  problem.      Suppose  now  that 
we  have  an  admissible,    normal   curve  which  furnishes  the  desired  maximum. 
On  it  the  Euler  equation  must  be  satisfied,    and  the  matrix    (I  *0   has   a 
zero  eigenvalue  with  a  corresponding  eigenvector   (clj    ;    we  will  choose 
c^>0.      Then  the   F  determined   by  the  Euler   equation  maximizes   the  pro- 
duct  K     BF,   as  a  function  of  0,   over  the  entire  interval   from  t  =  0  to 
t=T.      The  proof   is   as   follows: 

Suppose  that  9,(t)    is  the  argument  of  F  for  some   admissible  curve 
satisfying  the  Euler  equation  but  that  9-,(t)    does   not  maximize  the  pro- 
duct   K     BF(9)   on  some  sub interval    (t-pt2)   of   (0,T).      Then  there   is   some 
other  9,    say  92(t),    such  that   K*TBF(92)>  K^BFOj)   on  the   interval    (tx, 
t2).      Let  9   =  %2  on  tne   interval    (t^,ti    +  dt^)  ,    for  t^   +  dti<t2,   and 
let  6  9  =  fe|K1TB 1 1       +  e2K2TB-^F—     elsewhere  on  the   interval    (0,T)    . 
Then,    since  both  paths   are  admissible, 

(20)  dx  =  K1Tb[f(92)    -  FCexJ    dtx  +   I11e1   +   I12e2 
and 

(21)  0   =   K2Tb£f(92)    -   F(91)J   dtx   +   I21ex   +   I22e2« 

If  we  now  multiply  equation    (20)    by  c-i    and    (21)   by  c?   and   add,    remember- 
ing that   K*T  =  c1K1T  +  c   K2T,   we  get 
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(22)  Cldx  =   K^B    JF(62)    -    P(ex)J    d^ 

since!      lis   an  eigenvector  of   (ILJ)    corresponding  to   its    zero   eigen- 
value.     Since  the  right-hand   side  of   equation   (22)  and   the  constant  c, 
are  both  positive,   the  dx  of  equation   (22)    is  positive.      Furthermore   it 
is   possible  to   satisfy   equations    (20)    and    (21)    by   a  set  of   e-    which  will 
give  a  positive  dx  and   at  the  same   time  keep  dy  =  0.      To   see  this,  mul- 
tiply equation  (20)   by   c.    and  call   the   new  equation    (20  );    multiply 
equation  (21)    by  c2   and  call   the  new  equation  ( 21   )  •      Add   equation   (20   ) 

to   equation   (21*)    and   call   the  resulting  equation    (23).      Looking  now  at 

*T  IT  2T 

equations    (23)    and    (20),    remembering   that   K    '   =   CjK   '         c2K      ,    we  have 

(23)  cTdx  =   K*TB  [F(92)    -   F(91)J    dtx 
and 

(20)  0   =  K2TB  [f(92)    -   F(9X)J  dtj    +   I21ex   +   I22er 

*T    r  i  ?? 

But,    by  assumption,    K     B  LF(92)    -   F(9i)J>0,    and   I**   ^  0,   hence  we  can 

solve  these  two  equations   for  e,    and  e„.      Hence  it    is  possible  to   satis- 
fy equations    (20)    and    (21)    by  a  set  of  e-    which  will   give  a  variation 
which  will    in  turn  yield   a   larger  value  of  x(T).      Hence   if   9(t)    is   such 
that   the  product  K   TBF(9)    is   not   a  maximum  on  every  subinterval   of    (0,t), 
then  x(T)   will   not  be  a  maximum.      It    is  worthy  of  note  that  the   form  of 
the   conditions   for   an  extremum   is   independent  of   the   particular   problem. 

The  Weierstrass-Erdmann  corner  condition  is  an   immediate  consequence 
of   the  above  condition.      For  suppose  that   we  have  a  curve  whereon  the 
hypotheses  of  the  above  theorem  are  satisfied.      Suppose  further  that  tx 

is   a  point   at  which  the  control  variable  €Kt)    is  discontinuous.     Then 

it-,. 

(ciK   +  c9KZi)BF(0)      =  0.   Proof:  This  must  be  so:  otherwise  the 
1       2  ltx_ 

condition  just  found  would  not  be  satisfied  in  some  neighborhood  of  t. , 
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for   either  t>  t.    or  t<tn.      The  condition  that     K     BF(9)  |  =     K*TBF 

1  1  1— 

(0) L        is   commonly  known  as   the  Weierstrass-Erdmann  Corner  Condition. 

In  this  section  we  have  introduced  several  important  concepts  as 
they  pertain  to  a   linear  problem  discussed   in  the  following  iections. 

In  the  next  section  we  consider  a  combination  of  linear  systems, 
the  combination  being  non-linear. 
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III.      The  Weierstrass-Erdmann  corner  condition  for  a  more   general   problem. 

In  this  section  we  consider  a  problem  having  a  corner  at  a  variable 
time.  We  develop  the  Weierstrass-Erdmann  corner  condition  for  this  par- 
ticular  problem. 

Let  us   consider  the  differencial  equation 

(24)  X  =  A'X  +  b'f 

where  A*    is  A,    for  0<  t<  tj,    and  C,    for   t1<  t<T,    where  B*    is  B,    for 
0<t<t1,    and  D,    for  t1<t<T       X  is  the  matrix/*^;    F   is   the  matrix 
[    1\;    A,B,C,    and  D  are  2x2  matrices   of   functions  of   t  which   are  piece- 
wise   continuous   and  bounded   on  their   respective   intervals;    B  and   D  are 

non-singular.     T   is   fixed;    t-,    is  a  variable  to   be  determined.      The  con- 

o  2 

straint  on  F   is   as   before,    namely    (f-i)      +    (f~)      ■  1.      In  addition,   x 

and   y  must  be  continuous   at   t-,.      For  this  problem  the   admissible  curves 

are  the  allowable  curves   starting  at  X     and   ending  on  the  line  y(T)    =  yf. 

Note  that   this   problem   is  non-linear  on  the   interval    (0,T),    since 
A     and  B*    are  functions  of   t-.    as  well   as  of  t.      The  problem   is,    however, 
linear  on  each  of  the  two   sub- intervals    (0,t^)    and   (t,,T).      Hence  we  call 
this   a  semi-linear  problem. 

We  can  rewrite  equation    (24)   as 
|AX   +  BF,   0<t<  tj 


We  choose,  as  particular  solutions  to  K  =  -C  K,  tlta  solutions  K  (T)  =  LJ 
and  KZ(T)  =(,]•  For  t,<  t<T,  then,  each  K  is  a  function  of  T  and  of  t. 
Since  T  is  to  remain  fixed,   however,    we  shall   suppress   it  wherever   it 
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occurs  in  the  KL  and  shall  consider  the  KL  as  functions  of  t  only,  on  the 
interval  (t.,T) .   Let  us  define  a  set  of  particular  solutions  to  the 
equation  K  =  ATK  by  K^t-^)  =  Ki(t1+)  ,  for  i  =  1,2.   For  t  on  the  inter- 
val (0,tp,  then,  the  KL  are  functions  of  the  variables  t,  and  t  where, 
as  above,  we  suppress  m 

Before  continuing  further  let  us  make  one  substitution.  When  con- 
venient, we  shall  use  E  to  mean  B,  in  the  interval  (0,t..)  and  D,  in  the 
interval  (t,,T).   (Note  that  E  and  B*  are  the  same  matrix).  Using  this 

convention,  we  get 

T 

(27)  x(T)  =  K1T(t1,0)X(0)  +  \K1TEF  dt 

0 

and 

(28)  y(T)  =  K2T(t1,0)X(0)  +  \K2TEF  dt . 

0 

2       2 
Since  the  constraint  on  F  is  that  (f,)   +  (f2)   =  1 ,  we  can  replace  f, 

by  cos  9  and  £2   by  sin  0,  where  9  is  a  function  of  t.   Since  we  are  con- 
sidering both  0  and  t^  as  variable,  we  get,  from  equations  (27)  and  (28), 

(29)  dx(T)  =  ^x(T)   dt,  +  6x(T) 

d  tj     L 

and 

(30)  dy(T)  =  *y(T)   dt,  +  6y(T) 

0  1 1 

where  8x(T)  and  &y(T)  are  the  variations  in  x(T)  and  in  y(T)  due  to 
variations  in  the  control  variable  0. 

&x(T)  and  $y(T)  are  given  by  the  following  equations  for  sufficient- 
ly small  values  of  0: 

T 


■S 


VK1TE-||-  69  dt 


and 
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T 
(32)  6y(T)    =  Ck2TE-||-  66   dt 


IT    ^  F  2T    d  F 

Let   us  choose  the  special   variation  69   =  e^K     E* +  e-K     E  ■>  .        where 

en    and  eo   are  constants   which  remain  to  be   specified.      Then  equation   (31) 

becomes 

T 

(33)  6x(T)    =  ^(K1TE|^)(elK1TE|f  ♦   e^Eff-)    dt 

0 

with  a  similar  expression  for  6y(T).   Adopting  the  substitution 

T 

(34)  Iij  =  ^(KiTEH-)(KjTE-|L|_)  dt, 

0 
we  get 

(35)  6x(T)  =  e^11  +  e^12 

and 

21      22 

(36)  6y(T)  =  exI   +  e2I   . 

We  will  use  equations  (35)  and  (36)  later. 

The  variations  in  x(T)  and  in  y(T)  due  to  a  variation  in  t,  are 

more  complicated.   Rewriting  equations  (27)  and  (28)  to  show  the  way  in 

which  ti  enters,  we  get 

tX  T 

(37)  x(T)    =   K1T(t1,0)X(0)    +  \K1T(t1}t)BFdt    +   y1T(t)DF  dt 

o  tl 


and 


tl  T 

(38)  y(T)    =   K2T(t1,0)X(0)    +   \K2T(t  x  ,tj  BFdt    +  JK2T(t,  .IP  dt.- 

n  t- 


Hence 


1 
t 


IT  f        IT 

(39)  ^x(T)      dt,    =  |LK_(t1,0)x(0)dt1    +     \  ?K        (tn 

dti  L        £  t-.         x  1         J  dt,  1 

A  1 


1 


,t)BFdt,dt 


+  [(K1TBF)t         -(K1TDF)t     1    dt. 
tl_  Zl+J         1 


and   a   similar  expression  for    ^y(T) }    t^e  oniy  difference  being  that 

dti 
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IT  2T 

each  K       is   replaced  by  K     .      To   avoid  problems    in  notation   let   us  denote 

IT  IT 

by   6K       the  variation   in  K       caused   by  varying  t,    by  a  small   amount  dt    , 

IT        «fc  K 
i.e.    to   a  first-order   approximation  OK       =  -2-f; dt^.      With  this   nota- 

o    t^ 
tion,    equation    (39)    becomes 

(40)  |^(T)        dtx    =    6K1T(t1,0)X(0)    +  ^6K1T(t1,t)BF  dt 


[(K  TBF)ti_-(K  TDF)t  J    dtr 


Since   K     is   not   a   function  of   t,    on  the   interval    (t..  ,T)  ,  6K     =  0  on  that 
interval.      For  t  on  the   interval    (0,t-i),   6K     satisfies   the  adjoint  equa- 
tion,   and  hence  6K     =   -A^  K     on   (0,tv)-      Keeping  this    in  mind,    let   us 
consider  the   integral    in  equation    (40).      From  equation   (25),   BF  =  X  -  AX. 
With  this   substitution,    the   integral    in  equation  (40)    becomes 

tl  tl 

(41)  \6K1T(t1,t)X  dt    -     \bK      (t1,t)AX  dt. 

0  0 

Integrating  the  first  integral  by  parts,  we  get 


IT 


(42)  6K     (t,,t)X(t) 


-J(6K  * 


d^S^X  dt, 


But   6K     satisfies   the  adjoint   equation,   hence  the  integral    in   (42)    is 

IT  IT 

zero.      Hence  (42)    reduces  to  &K      (t1,t.)X(t1)    -6K     (tx,0)X(0).      But  the 

IT  ... 

variation  in  K       due  to  a  variation   in  t,  ,    evaluated  at   t-,,    is   equal   to 

LK1T(t1_)A  -  K1T(t1+)C/  dt,,    by  equation   (25).      Furthermore,   X(t)    is   con- 
tinuous,   so  that  X(t,    )    =  X(t1+).      Hence   (42)    is   equal   to 

LK1T(t1   )A  -  K1T(t1+)cJ   X(t1)    dtx   -  6K1T(t1,0)X(0),    and  equation    (40) 
reduces   to 

(43)  ^x(T)         dtj   =      K1T(tx_)    [AX(tx)    +  BF(t1_)J     dtx 
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-    K1T(t1+)    [cX(tx)    +    DF(t1+)J     dtj. 
But,   by  equation  (25),    this   says   that 

x(T)      dtx    =   [K1T(t1jX(t1_)    -    K1T(t1+)X(t1+)]      dtr 


(44) 

dx(T) 

Similarly, 

(45) 

5y(T) 

d    tn 

(T)  r    2T  .  2T  •  1 

^—  dtj-    =   [_K    i(t1_)X(t1_)    -    K      (t1+)X(t1+)J       dtr 

Hence  equations  (29)  and  (30)  give,  as  the  total  variations  in  x(T)  and 
y(T)  , 


T 
(46)     dx(T)  = 


and 

(47) 


(  K1TB||_6e  dt  -  (k1tx)  jtl+  dtx 

0  1_ 


T  t 

dy(T)    =  [   K2TE|2_  69  dt   -    (K2TX)    1    1+  dtj 

0  *Pu 


IT     2T  Ps  F 
Choose  the  special  variation  69  =  (e-,K  '  +  e~   )E?    .   Then,  using 

equations  (35)  and  (36),  we  get 


(48)     dx(T)  =  I11e1  +  I12e2  -  (K1TX) 


1  + 

dtl 
«1- 


(49)     dy(T)  =  I21ex  +  I22e2  -  (K2TX) 


Observe  that  if  t.  is  fixed  we  have  the  problem  discussed  earlier. 

Let  us  assume,  then,  that  admissible  curves  exist  and  that  our  prob- 
lem has  a  solution.   If  admissible  curves  exist,  it  is  possible  to  find 
e,,  e„  and  dt,  from  equations  (48)  and  (49)  and  hence  to  get  a  variation 
&9  such  that  the  curve  obtained  by  replacing  9  by  9  +  £9  is  admissible. 
Suppose  that  we  have  performed  these  calculations  and  have  obtained  an 

admissible  curve.   For  this  curve,  dy(T)  in  equation  (49)  is  zero.   Be- 

22 
fore  continuing  further,  let  us  assume  that  I/O.   This  assumption 
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is  enough  to  insure  normality  in  this  problem.   Having,  then,  an  admis- 
sible normal  curve,  we  can  use  the  results  obtained  after  equation  (19). 

We  showed  that  on  an  admissible  normal  curve  for  which  x  (T)  is  a 
maximum,  the  matrix  (I1-')  has  one  zero  eigenvalue  and  that  the  components 
c^  and  C2  of  the  eigenvector  can  be  chosen  so  that  c^  is  positive.  We 

also  showed  that  one  solution  to  the  adjoint  equation  which  gives  an  ad- 
it     1      9 
missible  extremal  is  K  =  c^K  +  c2K  .   Hence  let  us  choose  c-,  and  Cy   as 

components  of  an  eigenvector  of  the  matrix  (ILJ)  corresponding  to  the 

zero  eigenvalue,  with  c-^  0.   Multiplying  equations  (48)  and  (49)  by 

these  c^  and  adding,  remembering  that  the  KL  are  continuous  at  t  ,  gives 

*1+ 


=   K*T(tx)         X(t)  I  dtJ 


*!+ 


(50)  c1dx  = 

This   leads  to  the  following  theorem: 

THEOREM;     For  x(T)    to  be  a  maximum,    the  quantity  K*T(t1)      X(t) 
must  equal    zero. 

Proof:  If  the  above  quantity  is  positive,  any  dt,>  0  will  give  a  posi- 
tive dx(T)  and  hence  a  larger  x(T).  If  the  above  quantity  is  negative, 
any  dt^<0  will  give  a  positive  dx(T).  This  is  the  Weierstrass-Erdmann 
corner  condition  for  this   problem. 

The  next   section  uses  the  problem  we  have  been  studying  for  back- 
ground  and  considers   the  essentially  different  problem  obtained  by  intro- 
ducing  the  constraint  that  x(t^)   has   some  fixed  value.      We  develop  a 
numerical  routine   for  determining  the  solution  by  a  method  of  successive 
approximations. 
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IV.  A  numerical  routine  for  determining  the  maximum  x(T)  for  a  fixed 
value  of  x(t^) 

In  this  section  we  consider  the  problem  states  as  before  but  with 
the  different  constraint  that  x(t^)  has  some  fixed  value.   We  make  up  a 
curve  consisting  of  two  arcs,  each  of  which  is  an  extremal.   We  use  the 
method  of  variation  of  extremals  on  these  arcs  to  drive  the  resulting 
curve  to  admissibility  and  in  a  gradient  technique  to  determine  the  curve 
on  which  x(T)  is  a  maximum. 

Let  us  again  consider  equation  (24),  namely 

(24)     X  =  A'X  +  B'F 
with  the  same  conditions  as  in  the  beginning  of  Section  III  but  with  a 
different  constraint  on  t,,  namely  that  x(t^)  =  x, ,  where  x,  is  given. 

This  problem,  like  the  one  in  the  preceding  section,  is  semilinear. 
It  is  not  linear  on  the  entire  interval  (0,T),  since  A*  and  B'  are 
functions  of  t^  as  well  as  of  t;  it  is,  however,  linear  on  each  of  the 
two  sub-intervals  (0,t^)  and  (t-pT). 

The  differential  equation  we  have  been  working  with  is 
TAX  +  BF,  0<  t  <  t, 


(25)     X  = 


.CX  +  DF,  t,<  t<  T 


For  this  problem  the  admissible  curves  are  allowable  curves  satisfying 
the  above  differential  equation  and  such  that  X(<|)  =  XQ ,  x(t,)  =  x-,, 
and  y(T)  =  y~,  where  x  and  y_  are  given  constants.  The  adjoint  equation 
for  equation  (25)  is 


-CTK,  t,  <  t<  T. 

1  9  1 

We  choose  solutions   K     and   K*  to  the   adjoint  equation  such  that   K^T)    = 
/    )   and  K  'CD    =    I    )   and   such  that  B?  are  continuous   at  t   =  t,.      We  further 
define  the   2x2  matrix  X    as  before;    its   first   column  is  K1(t)    and   its 
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second  column  is  K    (t).      Note  that  X  is   a  function  of  the  variable  t   for 

t^<t<  T  and  of  both  the  variables   t^   and  t   for  0<t<t1. 

T 
Multiplying  equation   (25)   on  the   left  by  ^C   ,    integrating   from  t  =  0 

to  t  =  t,,   and  using  the  fact  that   7(  is   a  solution  to  the  adjoint  equa- 
tion leads   to 


(51)     7^T(t1,t)X(t)  |    =  ^(tpt)BFdt. 


0      0 
We  have  also,  from  the  last  section, 


f1  T 

,0)X(0)  +  \    K1T(trt)BFdt  +  )KlT 


(37)  x(T)  =  KiJ-(t1,0)X(0)  +  )    K  i(t1,t)BFdt  +  JK   (t)DFdt 

0  tl 

2T  C   2T  f  2T 

(38)  y(T)  =  K   (tx,0)X(0)  +  J  K   (t^OBFdt  +  \K   (t)DFdt. 

0  t^ 

The  differential  equation  is  linear  on  each  of  the  intervals. 
We  choose  F  by  using  the  method  of  variation  of  extremals  on  (0,t-,) 
and  each  of  them  (t,,T).   On  each  of  the  two  arcs,  we  have  the  following 

■k  ~k 

by  the  theorem  of  section  II.   If  C  is  an  admissible  arc,  if  K  is  the 

~k  1        9  1 

solution  to  the  adjoint  equation  defined  by  K  =  c,K  +  c„K  ,  where  K 
and  K2  are  the  solutions  to  the  adjoint  equation  defined  above,  if  c-,^  0, 

it         *     .  *T  * 

and  if ,  on  C  ,  F  maximizes  the  scalar  K  EF,  then  C  furnishes  a  max- 
imum x  at  the  end  of  the  arc,  relative  to  the  point  at  which  the  arc 
began.   If  we,  then,  choose  F  to  maximize  the  product  K  BF,  for  0<  t 
<tn,  and  K  DF,  for  t^  <t<T,  the  resulting  curve  will  be  made  up  of  an 
extremal  from  t  s  0  to  t  =  tj  and  an  extremal  from  t  =  t,  to  t  =  T. 

Once  we  have  an  admissible  curve  made  up  of  two  extremal  arcs,  we 
want  some  way  to  vary  these  curves  so  as  to  get  a  maximum  x(T) .  We  use 
a  gradient  technique,  explained  later,  to  choose  a  set  of  variations  in 
the  curve  parameters  to  get  a  larger  x(T) •   The  new  curve  may  not  be  ad- 
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missible  because  of   second-order  effects.      Hence,   we  again  drive  the 
curve  to   admissibility,    say,   by  the  method  of  variation  of   extremals 
mentioned   above.        We  continue  this   process   until   the  maximum  x(T)    is 
obtained. 

The  calculations  proceed  as  follows:  Since  c.  is  positive  on  each 
arc,  we  can  choose  c,  =  1.  Hence  K  =  K  +  cK  ,  where  c  is  a  constant. 
Unfortunately  when  a  condition  on  X   is  given  at   time  t^,    K     may  not   be 

continuous   at  t,.      Hence   let   us   suppose  that   c    is   not   the   same   in  both 

*  1  2 

arcs,    and   let   us   adopt   the   following   notation:      K     =   K     +   aK    ,    for  t  on 

(0,t^),    K*  =   K     +  bK^,    for  t  on   (t,,T),   where  a,   b,    and  t,    are   unknowns 

which  must   be  determined   so   that  the  curve  made   up  of   the   arcs    is   admis- 

JL  jl.  -fcT 

sible  and  yields  a  maximum  x(T).   Now  consider  F  .   F  maximized  K  BF 
on  the  interval  (0,t^);  on  that  interval  K  is  a  function  of  t,,t,  and 

JL. 

a.   On  the  interval  (t^,T),  K  is  a  function  of  b  and  t  only,  hence  F* 
is  a  function  of  b  and  t  on  the  second  interval. 

■k 

With   F  replaced  by   F  ,    equations    (51),    (37),    and   (38)    become 

(52)  -^T(t1,t)X(T)  =   $     7f  (tx,t)BF*dt 

V  0 

tl  T 

(53)  x(T)    =   K1T(t1,0)X(0)    f   U1T(t1,t)BF*dt   +   \K1T(t)DF*dt 


tl  T 

-2T,         rtN     ,_         C    2T,  „      *_  (    2T,    i       *. 

+    \K      (t1,t)BF  dt    +   \r 


(54)  y(T)    =   K      (t1,0)X(0)    +  jK      (tx,t)BF  dt    +  JK     (t)DF   dt 

0  H 

From  these   equations  we  want  to  devise   a  routine  that  will,    first,   give 

It    is  possible  that  we  might   obtain  an  x  (T)    from  a  set  of  curves 
that   are  admissible  under  one   admissibility  criterion  but   not   under 
another.      After  trying  various   admissibility  criteria,   the   author  de- 
cided upon  that  of  calling  a  curve  admissible  if    |x(t,)-x..|    +    |y(T)-yf\ 
$  10~\ 
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us  a  set  of  admissible  curves  and,  second,  find  an  admissible  curve  where- 
on x(T)  has  its  maximum  value.   We  want  also  to  find  a  condition  which 
will  indicate  that  no  admissible  curves  exist,  if  such  is  indeed  the  case. 
In  this  problem,  a  curve  can  be  characterized  by  a  set  of  values  for  t-,, 
a,  and  b.   Hence  we  want  a  routine  to  find  values  of  a,  b,  and  t.  which 
will  determine  an  admissible  curve  whereon  x(T)  is  a  maximum. 

Let  us  first  get  expressions  for  the  total  differentials  of  x(t1), 
of  x(T),  and  of  y(T) .   First,  of  x(t^).   Equation  (52)  can  be  written 

t"l 

(55)  ^CT(t1,t1)X(t1)  =  7CT(t1,0)X(0)  +J^(t1,t)BF*dt. 

Taking  differentials  of  both  sides  gives 

(56)  d^,T(t1,t1)X(t1)  +#T(t1,t1)dX(t1)  =  d-^T(t1,0)X(0)  + 

h 

^T(t1,ti)BF*(t1_)dt1  +  \    d  #T(t1,t)BF*dt  + 


tl 

^  ^T(t1,t)B6F*  dt. 


t. 

s 

0 

T 
To  get  d  /(  (t.,t  ),  note  that  the  first  t-,  denotes  the  end  of  the  inter- 


val (0,t-,)  and  that  the  second  is  the  value  that  the  running  variable  t 
assumed  at  the  end  of  the  interval.  The  differential  due  to  the  change 

IT1 

in  the  first  t,  is,  since /(  satisfies  the  adjoint  equation,  *ft   (t,,t  ) 
(A-C)dt^;  the  differential  due  to  the  change  in  the  second  t-^  is  -7£  (t^, 
t1)A  dt1.   Hence  d^T(t1,t1)  =  -  ^r(t1,t]L)Cdt1  and  d"^T(t  ,t1)X(t1)  = 
-7tT(t1,t1)GX(t1)dt1. 

f1    T 

\  d^(t1,t)] 


Consider  next  V  d  ^"(t, ,t)BF  dt.   By  the  argument  after  equation 


0 


(42),   this    integral   is   equal   to  H    (t-^t    )(A-C)X(t   )dt.,    -  d?ti(t    ,0)X(0). 
Substituting  the  above  results    into   equation   (56)    yields 
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(57)  ■^T(t1,t1)dX(t1)    =    ^T(t1,t1)AX(t1)dt1   ♦ 


fcl 


/CT(t1,t1)BF(t1_)dt1    +  j  /rT(tx,t)B6F*dt 


^^VV^l-^l    +S     ^T(t1,t)B6F*dt. 

0 

Calculating  the   integral    in  the  above  equation   is   somewhat  more  diffi- 

*  IT  2T 

cult.      F      is   chosen  to  maximize  the   produce   (K       +  aK     )BF.      Expanding 

IT  2T 

the  product    (K       +   aK     )B  yields   the  matrix 

(  [kU+ak12J  bn  ♦   [k21+ak22Jb21  tn^12Jb12.  t21.ak22;b22)    . 

Calling  the  first   element  h,    and   the  second  h~   and   remembering   both  that 
F     must  maximize  the  product    (h,    h2)F  and   that  F  can  be  written  as 

(cos  9V   we  see  that  cos  0  =  h-t/   [(hx)2   +   (h2)2J  *  and   that   sin  6   =  h2/ 
sin  0/  . 

r         ?  91k  *        /-sin  8\  _1 

[(hx)      +    (h2)zJ    2.      Hence  6F     =   (  J  60,   where  0  =   tan        (h  /hj). 

Substituting   in  the  expressions   for   sin  9   and   for  cos   0   and   finding  60 
in  terms  of  6  h^   and   6h2  yields 


-h2\        h16h2   -  h26hx 
V    [(hx)2   +(h2)2J3 


6F*   =    , 


6h,  and  6 h2  are  each  the  result  of  variations  both  in  t  and  in  a.   Con- 
sider first  the  part  of  the  variation  (h  6h  -h„6h  )  due  to  a  variation 


in  a. 


*h9  3h. 


(58)    hjsp!  -  h*.  =    [(k11+ak12)bu   +    (k21+ak22)b2]]  [k12b12   +  ^^^ 

-[(kU+ak12)b12   ♦   (k21+ak22)b22][k12bu+k22b21J 
=   a  [(k12b11+k22b21)(k12b12+k22b22)-(k12b12+k22b22)(k12bn 
+k22b21)]      +[(k11b11+k21b21)(k12b12+k22b22)    -    (kUb12 
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21  12  22        "I 

+k     b22)(k     bll    +  fc      b2J 


=  k11k12(b11b12-b12bu)    +  k     k      (b21b22-b22bn) 

11   22  21    12 

+k     k      (bnb22~b12b21')    +  k     k      (b21b12~b22blP 

=    (k11k22-k21k12)(bnb22-b12b21) 

;   i-|B|, 

where   1 7\. |    denotes   the  determinant  of    /C  and   similarly  for    IB    .      Hence 

the  variation  due  to   a  variation   in  a   is  \K.\    |  B|    da. 

The  part  of  the  variation   (h^6h2-h26h-^)    due  to   a  variation   in  ti    is 

^h2  ^hl  T    11        12  21        22  -l 

(59)         <hiT£-  -   h2lt^    dtl   =     L(k      +ak      )bll  +  (k      +ak      )b2lJ 

•[(6k11   +  a6k12)b12+   (6k21+a6k22)b22]- 

jCk11   +   ak12)b12   +    (k21+ak22)b22J  £(6kU 

+a6k12)bu   +    (6k21   +  a6k22)b21] 


=  B 


k  +  ak     ~6k  +   a6k 

21    22    ,  21    22 
kZi+  ak     ^k  +  a6k 


For  convenience,  we  shall  refer  to  this  variation  as  VAR(l)  henceforth. 

i  i  1/  T  T 

Note  that  each  of  the  6k  J  can  be  calculated,  since  6/(,(t,,t,)  =  (A  -C  ) 

jf^(t1,t1)dt1  and  since  b"K,=   -AT6^on  (0,^). 

It  is  possible  to  simplify,  within  the  integral,  the  terms  /£  (t-,,t)B 

(21,   as   follows: 


/-h 
X  (tl.t)B 


k11   k21\  /bxl   b12\    /-    (k11+ak12)b12+(k21+ak22)b22 
k12  k22/  lb21   b22 J  V     (ki:L+ak12)b1:L+(k21+ak22)b21 


■("JHW 
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Using  the   above  substitutions,    equation  (57)    simplifies   to 
(60)  -^T(t1,t1)dX(t1)    =    ■^T(t1,t1)X(tl_)    + 


h 


a W|    If1^-^^     at. 

1/  l(h1)2+(h2)2J  3/2 


"fCT  is   a  non-singular   square  matrix,    hence   (*K   )~     exists,    and 


(61)  dX(tx)    =  X(t1_)    dt1    +    \fC    (t    ,tx)J 


-1 


[  (~%W        W  Id*.  -  va»(1)  at,      dt> 

where  dX(t,)    =  (  \  )    •      For  future  reference,    let   us   rewrite  dx(t.) 

1  \dy(t1)y  *■ 


as 


(62)  dx(t1)   =  a11dt1   +  a^da, 

£x(tx)  >r(tX) 

where  a.,    =  — r and   a , .   =  — .      Thus   we  have  derived   an  ex- 

11  ©>    t1  12  ^    a 

pression  for  dx(t,)    in  terms  of  variables   which  we  can  calculate. 

We  want    next   to  derive  an  explicit   expression   for  dx(T).      We  can 

rewrite  equation   (53)    as 

T 

(63)  x(T)    =    K1T(t1,t1)X(t1)    +  J    K1T(t)DF*dt, 


*1 


whence 


(64)  dx(T)    =   dK1T(t1,t1)X(t1)    +   K      (t1,t1)dX(t1) 

T 
IT  *  /"IT  * 

-   K      (t)DF    (t      )dt      +   \    K      (t)D6F  dt . 

tl 

IT  IT 

By  the  discussion  preceding  equation  (57),  dK   (t-,,t,)  =  -K   (t  ,t,)Cdt,, 

and  dK1T(t1,t1)X(t1)  =  -K   (t1,t.)CX(t1)dt1.   If  we  write  the  differen- 
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tial   for  dy(t.)    corresponding  to  equation   (62)    as  dy(t,)    =  a      dt      + 
ou^da,    the   product   K1T(t1,t1)dX(t1)    becomes 

11  12  11  12 

(kiJau  +  ki/a51)dt1  +    (k     a12  +  kiZa52)da. 

•k 

By  an  argument   similar  to   that   preceding  equation   (60),    noting  that  F 

is  dependent  only  on  b   in  the   interval    (t^,T),   we  get 
T  T 

(65)  \    K1T(t)D6F*dt   =  j    (-b)       IbI  2   (jtPdt1db 

H  fci     ^2,(^)233/2 

Combining  the  above  results  yields 

(66)  dx(T)    =  -K1T(t    ,t   )  fcx(t1)+DF*(ti    )]dt.    +  (kUa     +k12a     )dt. 

11""        L  1+ •»      1  11  ->i        i 


2dt  db. 


+  (k11a19+k12a£  Jda+  \(-b)      [b|  2  \Kl 

62   I  [(h1)5i(h;,2j3/2 


Combining  terms  gives  us   an  equation  corresponding   to   equation    (62), 
namely 

(67)  dx(T)    =  a21dt1   +a22da  +    a23db»   where01^  =    ^x(T)      , 

a„   =    |£fi3      ,    and  a32  =   SjS^  . 
22         }    a  d    b 

dy(T)    can  be  derived  much   as   we  derived  dx(T);    the  equation  corresponding 
to   equation    (66)    is 

(68)  dy(T)    =   -K2T(t1,ti)[cX(t1)    +   DF*(t1+)Jdt 

+   (k21a11+k22a^1)dt1  +   (k21a12+k22ag>2)da 

T 
+  C    lDf2lft|2   dt  db, 


H 


[(h1)2+(h2)2j3/2 


which  we  abbreviate  as 


(69)  dy(T)    =a31dt1    +  a     da   +  ot33db. 
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We  thus  have  three  equations   to  use   to  determine  corrections    in  t,,    a, 
and  b  which  will   lead  to   a  maximum  x(T),    namely 
(62)  dx(tx)    =  <xudt      +  a12da 

(67)  dx(T)      =a21dt1   +  a22da   +a23db 

(69)  dy(T)      =a31dt1   +  a32da   +  a^db 

The  numerical    procedure   used   for   solving   for   a  maximum  x(T)    is   the 
following:     Choose  values   for   t,,    a,    and  b.      With  these  values,    calcu- 
late x(tx),   x(T)  ,    and  y(T)    and   the  o(  .  .,    starting   from  X(0)    =  X   .      Prob- 
ably the  curve  so  determined  will  not  be  admissible,    i.e.    x(t   )    will   not 
be  x,    and  y(T)    will    not  be  yf.      To  obtain   a  curve  which  will  be  admissi- 
ble,   solve  equations   (62)    and    (69)    for  two  of   the  variables   dt    ,    da,    or 
db,    setting   the  third   equal    to    zero.      It    is   possible  to   solve  equations 
(62)    and    (69)    for  dt..  ,   da,    and  db  only   if  the  rank  of   the  matrix 

(70)  N     ""        °    \ 
\a31     a32      a33' 

is  two,  hence  this  is  a  criterion  for  deciding  whether  any  admissible 

curves  exist.   If  the  rank  of  (70)  is  two,  replace  a,  b,  and  t  by  a  +  da, 
b  +  db,  and  t-.  +  dt   and  repeat  the  calculations.   Continue  this  process 
until  either  the  rank  of  (70)  becomes  one  or  until  some  admissibility 
criterion  is  met.   Note:   In  correcting  two  variables,  it  is  essential 
that  the  value  of  third  variable  be  chosen  close  to  the  value  it  will 
have  on  an  extremal;  otherwise  the  correction  routine  may  not  converge 
to  an  admissible  path,  even  when  such  a  path  exists.   The  admissibility 
criterion  used  by  the  author  was  Ix-^  -  x(t-,)|   +   lyf  -  y(T)|  $   10~  ;  for 
the  matrices  used  this  criterion  was  normally  met  in  less  than  ten  iter- 
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at  ions  when  the  determinants 


Pll  a12 


|a31  a32 


and 


all     0 
*31    ^33 


went  to    zero,    this 


was   also  obvious  by  the  tenth   iteration. 

When  admissibility  had   been  achieved,    equations    (62),    (67),    and 
(69)   took  on  the  form 

(71)  0  =  a11dt1   +  a12da 

(72)  dx(T)    =  a21dt     +  a22da  +  <x23db 

(73)  0  =  a~,dt     +  ouoda  +  a33db. 

These  equations  were  solved   for  dx(T)    unless  the  rank  of 

ail   «12  \ 


a21   a22   a23 


/an   0^2   0    V 

was   equal    to   the  rank  of  \   .      If   the  rank  of 

v  a       a       a 
\  31      32      33/ 


^31    a32   a33/ 

the  latter  matrix  was   two,    this  condition  was   that   the  determinant  of  the 
first  matrix  must  be  zero.      We  therefore  want   some  way  to  determine  a  new 
set  of  corrections  to  drive  the  determinant   toward   zero.      The  method 
adopted  was   to   choose  dt,,   da,    and  db   proportional   to   the  components  of 
grad  x(t..)   X  grad  y(T)  ,   where   "X"   indicates   the  vector  cross-product   and 
where  the  proportionality  constant   is  some  fixed  number  e    multiplied  by 
the  determinant    ja.  .  I   .      The  old  values   for  t,,a,   and  b  were  replaced  by 
the  corrected  values    and   a  new  admissible  curve  was   found   using  the  meth- 
od outlined   above.      For  this  curve   the  determinant   |ttj.:\    and  the  correc- 
tions  da,    db,   and  dt,    were  again  calculated   together  with  x  (T) .      This 
process  of   finding   an  admissible  curve,    then  calculating  a  set  of   cor- 
rections  to    increase   x(T) ,   then  finding  a  new  admissible  curve,   then 
finding  a  new  set  of  corrections  to   increase  x(T)   was   continued  until 
the  determinant     ^x.  .\    changed   sign.      When  the  determinant  changed   sign, 
however,    the  above  procedure   no    longer  worked   -  the  corrections  became 
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larger  rather  than  smaller,  and  the  values  of  x(T)  got  smaller  instead  of 
larger.   A  procedure  that  worked  when  the  determinant  changed  sign  was 
the  following:  A  new  proportionality  constant  was  chosen  equal  to  one 
fifth  of  the  sum  of  the  positive  and  negative  determinants;  a  new  set  of 
corrections  da,  db,  and  dt-,  was  then  calculated.   The  admissible  curve 
calculated  after  finding  this  new  set  of  corrections  invariably  yielded 
a  larger  x(T)  than  had  been  obtained  before.   This  method  of  successive 
approximations  seems  to  be  a  nearly  foolproof  method  for  determining  the 
maximum  x(T) .   The  procedure  was  stopped  when  either  the  absolute  value 
of  the  determinant  was  less  than  10~5  or  it  became  obvious  that  the  con- 
vergence was  too  slow  for  the  amount  of  time  available  on  the  computer. 
If  convergence  was  too  slow,  it  could  usually  be  improved  by  increasing 
the  number  £  referred  to  above. 

Another  successive  approximation  procedure  used  was  to  compare  the 
x(T)  resulting  from  replacing  t^,a,  and  b  by  their  corrected  values  with 
the  previous  x(T).   If  the  new  x(T)  were  smaller,  £  was  reduced  by  a  fac- 
tor of  ten  and  a  new  set  of  corrections  and  hence  a  new  x(T)  computed. 
This  process  was  continued  until  either  the  sum  of  the  absolute  values 
of  the  corrections  was  less  than  10"* 5  or  the  new  x(T)  was  larger  than 
the  old  x(T) ,  in  which  case  the  new  x(T)  became  the  value  compared  with. 
This  process  gives  a  monotonically  increasing  sequence  of  values  of  x(T); 
the  speed  of  computation  was  comparable  to  that  of  the  other  method. 

It  was  necessary  to  insure  that  the  corrections  from  the  routine 
were  not  so  large  as  to  make  the  linearity  approximations  in  the  correc- 
tion integrals  invalid.   The  author  used  the  following  procedure:   In 
the  correction  to  admissibility,  if  either  da  or  db  exceeded  .3  in  abso- 
lute value,  or  if  dt-^  exceeded  .5  in  absolute  value,  he  divided  all  the 
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corrections   in  half.      For  the  problems   worked,    this   criterion  was   ade- 
quate although   in  several    cases    it   slowed  convergence  quite  a  bit.      (For 
example,    where  the  initial  value  given  for  a  was    1.    and  where  the  final 
value  on  the   admissible  curve  was   12.8.)      Once  admissibility  was   achieved, 
it  was   necessary  to  use   some   linearity  criterion  on  the  second   set  of 
corrections.      The  procedure  was   to   divide   the  corrections   in  half   if  the 
sum  of  the   absolute  values  of  the  corrections   exceeded   .5. 

With  the  above  correction  procedures,    computations   were  made  on 
several   sets  of  matrices   with  various   boundary  conditions.      The  only  case 
where  it  was   possible  to  verify  the  results  by  comparison  with  a   physical 
situation  where  A  and  C  were  set  equal   to    zero   and  where  B  and  D  were 

2x2   scalar  matrices.      This   set  of  matrices  gives   the  differential   equa- 

("aF,    0  <t<tx  ?? 

tion  X  =  1  t  +  *  r  '      ^he  constrai-nt  on  F  was   that    (f-^)      +    (f~)      =   1; 

f-i    and   f2  were  consequently  chosen  as   cos   a  and   sin  a,    for  0<  t<   t,  ,    and 
as   cos    B  and   sin  B,    for   t,  <  t  <  T.      Solutions   to   the  differential   equation 
as   set  up  describe  the  path   taken  by  a   light   ray  going  from  one   isotropic 
medium  into   another,   with  a  and    B  being  the  angles   between  the   light   ray 
and  the  normals   to   the   plane  of  discontinuity.      As   such,    these   solutions 
should  obey  Snell's  Law,    namely  that   the  ratio   of  the  velocities  on  oppo- 
site sides  of  the  discontinuity  plae   should  be  the  same   as   the  ratio  of 
the  sines  of  the  angles   the  ray  makes  with  the  normal   to   that  plane.      The 
cases   tried  were  for  X     =   (n)>    x,    =   1.0,   y~  =   5.0,    a  =   1.0   and  b  =   2.0 
and,    in  the  second   case,    a  =   2.0   and  b   =   1.0.      For  the   first   case,    the 
ratio   between  the   sines   should   be    .5.    a    was   found  to   be  approximately 
.403378   radians   and  A  approximately    .902758   radians.      The  ratio   of  the 
sines  of  these  angles  was    .500,   which  verified  the  accuracy  of  the  rou- 
tine.     For  a  =   2.0   and  b  =   1.0,   the  angles   were  reversed,    and  again  the 
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accuracy  of  the  routine  was  verified. 

A.O   ,l\      A.O   .2\ 
Another  case  tried  was  where  A=  I      nJ'B=\   oil/' 

/3  0\  /2.0  0.  \ 

C  =  [     J,  and  D  =  /     ?  oi '  witn  xl  =  1,0»  T  =  **#0»  and  where  initial 

values  for  t,,  tana,  and  tan  0  were  .73,  4.0,  and  .7,  respectively. 
The  routine  converged  to  tj_  =  .8935,  tan  a  =  1.4925,  and  tan  p  =  .4967, 
with  the  maximum  x(T)  being  6.4622.   This  convergence  took  fifteen  min- 
utes and  fifty-one  seconds  on  a  CDC  1604  computer.   The  author  felt  that 
convergence  could  be  improved  if  better  initial  estimates  were  made  of 
t-,  ,  tan  a,  and  tan  0. 
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APPENDIX 


O 


These  constants  are  given: 

x(0).  y(0),  XMID.  YFINAL,  MATRICES  A, 
B.  C,  and  D,   EP8FAC,  Xll,  ATCi2,  ATD11, 
TONE,  CGUES8,  DGUES8,  CAPT,  ERROR. 
EPS1 


SET  STEP  SIZE 


TSTEP(l)  =  (CAPT  -  TONE)/50 
TSTEP(2)  =  TONE/50 


I 


INTEGRATION  FROM^  ■  CAPT  to  t_°_TON-E 
In: /f (CAPT)    =  (*  °) 

Integrate:  K  ■   -CT/T 
Out:  W(t~)  ~ 


I 


Calculate:  6*-^)    =   (AT  -  CT)A"  (tj) 


INTEGRATION  FROM  t 

=  TONE  to  t  =  0 

In:  A*^). 

O/^tj) 

Integrate: 

K  =    -ATA, 

6/r 

=   -AT6K  ' 

Out:  K  (0).    6/T(0) 


INTEGRATION  FROM   t  =  0   to   t  =  TONE 

to!(y(0))  '  *(0)'   6/r(0>'       Note:  *  =  (K*:K^  :     6/r  *   (oK^flK2) 


-ATA\    6* 


T 
-A*6K 


Integrate:    K 

Calculate:     (aCON(1)     ACON(2))     =    (k1T   +  CGUESS   x   K27)   B 

SRT   =    MACON(l))2     +     ^ACON(2A2^ 

F     =     (aCON(2))    *  SRT 


Integrate: 


AX    +   BF 


d    (cORR(l))    "  (-CGUESS   x    |b(2    x  JA'|  2     4  (SRT)3)  dt 
d    (cORR(2))    =  (|B|2   x    \k\2       +(SRT)3)dt 


Calculate:      DINTFAC   =     |B| 


K11    +   CGUESS   x   K12 
K21    +   CGUESS   x   K22 


NOTES 


I. FLOW  CHART  PAGE  NUMBERS  ARE 
DESIGNATED   BY  THE   LETTER   F.  e.g., 
FG    IS   PAGE  S  OF  THE   FLOW  CHART. 

2.C0NNECT033  ARE    REFERENCED    BY 
PAGE  NUMBERS,  t.g., (7^  F2  DESIGNATES 
CONNECTOR  3  IS  ON  v-^  PAGE   F2.  J 


6KU    ♦  CGUESS  x   6K12 
6K21    ♦  CGUESS  x   6K22 


Integrate:       d     (CORR(3))    =     -CGUESS   x    |b|2    x    |/f|   x   DINTFAC    +   (SRT)3  dt 
d     (C0RR(4))    =     |B|2   x    \K\    x   DINTFAC    *  (SRT)3  dt 


x(t,) 


1 


I 


Calculate   CKI   =     (^i^))'1 


I 


CHANGES  IN  X(tj)  AND  Y  DUE 
TO  CHANGES  IN  CGUESS 


XONE(2)    =   (1  0)  CKI 


/CORR(l)\ 
VcORR(2)/ 


YONE(2)    =   CORR(2) 


I 


VARIATION  IN  X(TONE)  DUE  TO  A  VARIATION1 
IN  TONE 

XONE(l)    =   (1  0)  CKI  (cJSrJJ))    +   XDOT1 


I 


'PARTIAL  VARIATION  IN  Y  DUE> 
TO  A  VARIATION  IN  TONE 


YONE(3)     =     CORR(4) 


0R 


/XDOTl\ 


ly(t)).     *  (*i-)»    X(tl-)     =    VydOTi)'    CORRW.    CORR(2),    CORR(3) ,    CORR(4) 


I 


Fl 
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£_ 


INTEGRATION  FROM  t  =  TONE   to   t  =  CAPT 
In:      tfCtjJ.     X(tx) 

I  T 

Integrate:      K  =   -G  K 

Calculato:    (ACON(3)   ACON(4))   -  (k1T   +  DQUE8S  x   K27)  D 

SRT  «=HACON(3))2    +  (aCON(4))2J 

_      /aCON<3)\  .  SRT 
F-      ^ACON(4);t5RT 

Integrate:     X   =  CX   +   DF 

On  first  step  only,  set  (^cm)    =  cx   +  DF 
Integrate:  /  ■  2  j  a\ 

D(CORA)    =  (  -DGUESS   x    \o\     x  l/Cl  "     +    (SRTT )  dt 

EKCORB)    -[   |D|2   x    |A"|2    +    (SRT)3)dt 


Out:  x(T),  y(T),  CORA,  CORB,  XDOT2,  YDOT2,  K  (T) 


XBNO  =  x(T) 
YBND  =  y(T) 


VARIATIONS  IN  X(CAPT)  AND  Y(CAPT)  DUE 

TO  A  VARIATION  IN  TONE 

XONE(5)   =   CORR<3)    ♦  KX\)  ("J™   I   ggg) 
YONE(l)  =   YONE(3)    ♦  K*^)  (jjgg    I   *jgg) 


Changes  in  X(CAPT)  and  Y(CAPT) 
due  to  change  in  DGUESS  are 
L  CORB  and  CORA  / 


XDIFF   =  XMID   -   XINMDJ 
YDIFF   =   YFINAL   -   YBND 


I 


£ 


/*   Path  admissible  ? 
\~  JXDIFF |  + 


|DELTATJ    >  2? 
TNo- 


Yes 


Set  DELTAT   =   DELTAT/2 


Print  DELTAC,  DELTAT, 
DELTAD 


O" 


TTWO   -   TONE    +   DELTAT 


c 


J 


Yes 


TTWO  so? 


No 


y± 


TONE    =   TONE/2 


± 


CGUESS   =   CGUESS    ♦  pELTAT  X  DELTAC 
DGUESS   =   DGUESS    +  r^ELTAT  x  DELTAD! 


0" 


c 


TTWO    a  CAPT  ? 


"> 


No 


Yes 


V 


TONE    =   TTWO 

CGUESS  =  CGUESS  +  DELTAC 

DGUESS  =  DGUESS  +  DELTAD 


1 


TONE 


TONE  +  CAPT 


I 


6 


F3 


nrnvaa  .  rrnir«ss   *  CAPT  -  TONE  .  DELTAC 

CGUESS   -  CGUESS   +  -  x  D£LTAT 

nriipco  -    rv-nrts   4.  CAPT  -  TONE  .  DELTAD 

2  x  DELTAT 


_0F3 


|YDIFF|    s  ERROR 


Jp-O 


DETER 


XONE(l)      YONE(l) 
XONE(2)     CORR(2) 


(   8JF4 


39  (F2) 


(cat 


CAPT   -   TONE    s   .0001  ? 


G 


No 


TONE    £  .0001  ? 


I 


No 


Print: 


CGUESS,  TONE 


Ofi 


Yes 


>* 


Print:    TONE 
TOO  LARGE 


Yes 


XONE(2)  0  XONE(l) 

CORR(2)       CORA      YONE(l) 
CORR(l)       CORB      XONE(5) 


Print:    TONE 
TOO  SMALL 


9 


Calculate   Un"1  (CGUESS)  and  Un"1  (DGUESS) 
ATANC    -   ATANF  (CGUESS) 
ATAND   =   ATANF  (DGUESS) 


Print  ATANC,  ATAND 


DELTC    -   ATANC    -   ATC11 
DELTD   «   ATAND   -   ATD11 


a 


i 


DELTcl      t     IdELTdI     <  10"6  ? 


No 


ATC11.  »   ATANC 
ATD11    -   ATAND 


o« 


Yes, 


Print:  WE  ARE  DONE 


EPS   =   EPSFAC   x    DTMT 


^* 

(x(T) 

Ve 

>    xn  ?          J      r 

ANo 

x(T)     is  new    Xll 

^ 

T 

' 

Calculate    DELT1T,  DELTIC,  DELT1D   so  as  to 

improve  x(T)  for  admissible  paths, 

EPS 

=     EPS     x     .1 

DELTIC    =    EPS    x  I       °                      XONE(l)| 
DELTIC    =    EPS    x  |  CQRA                 YONE(l)| 

«»».»., |~g      g-gj 

»»»••«■«  "IS!       coJa    1 

\ 

CGUESS    = 
DGUESS   = 
TONE 

CGUESS    -    DELTIC 
DGUESS    -    DELT1D 
TONE    -   DELTAT 

DELT1C(      +     |DELT1D| 
+  IdeltitI    <    .5  ? 


Yes 


,No 


Print    DELTIC,    DELT1D,    DELT1T,    EPS 


Divide  all  corrections 
by  two 


►QF3 


CGUESS 

= 

CGUESS 

+   DELTIC 

DGUESS 

* 

DGUESS 

+   DELT1D 

TONS 

= 

TONE    ♦ 

SELTlT 

-o 


F3 


40  (F3) 


MdeterIs    10~5    7      ) 


No 


Solve  for  DELTAD  and  DELTAT 


/xdiff\       /xone(i)        o    \/deltat\ 

^YDIFFJ    =    ^YONE(l)       CORA  /  \DELTAiy 


©■ 


-> 


Solve  for  DELTAC  and  DELTAD 

/XDIFF^    „    /XONE(2)  0 

\YDlTFj         ^CORR(2)       CORA 


)(: 


DELTAC \ 

DELTAD ) 


>  (JDELTACl 


G> 


■v  No 
<3?        -H 
Ye8   £ 


Divide  all  correc- 
tions by  two 


10  )F4 


(iDELTADl     <    .3     ?  )  Ye8  » ^lxONE(2)   x   CORR(2)|  *  lo"5    7  \^-»>(|DELTAT|        <  .5     ?    Y"     frfP) 

i  ■  i  Yes 


F2 


Divide  all  corrections  by 
two 


6 


F4 


Print:    NO  ADMISSIBLE 
PATHS  EXIST 


Divide  all  corrections  by 
two 


Cll\  F4 
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APPENDIX   II 


PROGRAM    DISCON  00f 

C  THIS    PROGRAM    COMPUTES    MAX    X(CAPT)     FOR    A    GIVEN    Y(CAPT)    AND    A    GIVEN  001 

C  X(T1)»WHERE    XDOT»AX+BF,FOR    T    BETWEEN    0    AND    T1»AND    CX+DF,FOR    T 

C  BETWEEN    Tl    AND    CAPT.    CAPT     IS    GIVEN, BUT    Tl     IS    NOT. 

DIMENSION    A(2t2)»B(2»2)»C<2»2),D(2,2)»BK(2»2).F(4) ,TSTEP(2 ) tE(2,2» 

1  2) >ELV( 14) »DE<16)  .         AAK ( 2  *  2 ) »DP ( 16 ) » ACONJ 4) . AKTB( 2 » 2 ) ,CORR  OOl 

2  (4),CKI(2.2),FLES(2) »XONE(5) »Y0NE(4) ,BKA(2.2) , FMOR ( 2 ) »DELK(2»2)  OOi 
3»DACON(4) ,DFE(4) ,AK(5»16)  00» 

READ  1503,EPSFAC  00» 

1503  FORMAT ( 1F20.10)  00» 

READ  501»( ( A( I,J) »J»1,2) tl*lt2) ♦< <B( I,J) ,J*1,2) ,1=1.2)  00 

501  FORMAT  (8F10.1)  00 
READ  501, ( (C( I »J) tJ=l»2)»I=l»2)  ♦(  <DU »J) ,J=1»2) » 1=1,2)  00 
PRINT  551  GO! 

551  FORMAT  ( 14X , 1HA,29X, 1HB » 29X , lHC »29X , 1HD )  00 
PRINT  5  52»UA(I,J).J  =  lt2)»(B(I,J),J  =  l,2)»(C(I,J),J«l,2)»(D(I,J),J  =  00 

•  +  1,2)  ♦  I  =  1 »  2 )  00 

552  F0RMAT(4( 8X ♦ F4. 1 ,6X » F4 . 1 ♦ 8X ) )  00 
READ  502, TONE»X0,Y0,YFINAL,XMTD»CGUFSS, ERROR, CAPT  00 

502  FORMAT  (8F10.8)  00 
READ  503,DGUESS  OOi 

503  FORMAKF10.8)  00. | 
PRINT  553, TONE, X0,Y0,YFINA|_,XMID,CGUESS, ERROR, CAPT  00; 

553  F0RMAT(4X,12HT0NE  GUESS=  »F5.2»  9H ( X0» YO ) = ( , F4. 1 , 1H , ,F4. 1 ,  9H)YFIN  00; 
+AL=  »F4.1»7H  XMID=  »F5.2»9H  CGUESS«  .F4.1.8H  ERROR*  »F7.6»6H  CAPT*  00; 
+,F5.2)  oo; 

C      INTEGRATE  ADJOINT  BACKWARDS  TO  GET  K(0)  AND  KlNVERSE  AT  Tl.  00; 

F(1)=0.  oo; 

F(2)=0.5  00. 

F(3)=0.5  oo; 

F(4)=l.  00' 

ATC11=5.  00' 

ATD11=5.  00' 
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Xll=l.  0033 

EPS2=-.01  0034 

TSTEP(2)=TONE/50.  0035 

DO  1003  I =1 »3  0036 

YONEU}=0.  0037 

XONE(I)=0.  0038 

TSTEP( 1 )=(CAPT-TONE)/50.  00^9 

TSTEP( 1)=-TSTEP(1)  0040 

TSTEP(2)=-TSTEP(2)  0041 

BK(1»1)=1.  0042 

BK(2,1)=0.  0043 

BK(1»2)=0.  0044 

BK(2,2)=1.  0045 

T=CAPT  0046 

DO  101  1=1,2  >                                                                       0047 

DO  101  J=l,2  0048 

E(2,I,J)=A( I ,J)  0049 

101  E( 1»I,J)=C( I »J)  0050 

ELV(1)=BK(1.1)  0051 

ELV(2)=BK(1»2)  0052 

ELV(3)=BK(2»1)  0053 

ELV(4)=BK(2,2)  0054 

DO  1004  IMA=5t8  0055 

ELV(IMA)=0.  0056 

IB=1  0057 

DO  102  IE=1»50  0058 

DO  103  IC=1,4  -                               0059 

DO  104  ID=1,4  0060 

104  DE( ID)=ELV( ID)+F( IC > *AK ( IC-1 , ID)  0061 

DP(1U-E( IB»1»1)*DE( 1)-E( IB,2,1)*DE(3)  0062 

DP(3)=-E( IB,1,2)*DE( 1)-E( I B,2 ,2 ) *DE< 3 )  0063 

DP(2>=-E( I8,1,1)*DE(2)-E( IB,2,1)*DE(4)  0064 

DP(4)=-E( IB,1»2)*DE(2)-E( IB, 2 ,2 ) *DE(4 )  0065 

DO  103  ID=1,4  0066 

103  AK( IC,ID)=TSTEP(IB)*DP(ID)  0067 

DO  105  ID*lt4  0068 
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105 
102 


1002 


904 


903 

905 
112 


ELV( ID) 

T=T+TST 

BK( ltl J 

BK( 1.2) 

BK(2»1) 

B<(2.2) 

IB  =  2 

ELV(5)= 

ELV(6)= 

ELV(7)= 

ELV(8)= 

DO  112 

DO  903 

DO  904 

DE( ID)  = 

DP(1)=- 

DP(3)=- 

DP<2)=- 

DP(4)=- 

DP(5)=- 

DP(7)=- 

DP(6)=- 

DP(8)=- 

DO  903 

AK( IC.  I 

DO  905 

ELV( ID) 

T=T+TST 

BK(l.l) 

BK(1.2) 

BK(2,1) 

BK(2,2) 

DELKU, 

DELKd. 

DELK<2. 

DELM2, 


=ELV( ID)+(AK< l.ID)+2.*AK(2»ID)+2.*AK(3»ID)+AK(4, ID) )/6. 

EP( IB) 

=  ELV(1  ) 

=ELV(2) 

=ELV(3) 

■ELV(4) 


(A(  1 
(A(  1 
(A(  1 
(A(l 
IA  =  1 
IC=1 
ID=1 
ELV( 
E(  IB 
E(  IB 
E(  IB 
E(  IB 
E(  IB 
E(  IB 
E(  IB 
E(  IB 
ID  =  1 
D)=T 
ID  =  1 
=  ELV 
EP(  I 
=  ELV 
=  ELV 
=  ELV 
=  ELV 
1)=E 
2)«E 
1)»E 
2)*E 


»1)-C( 
♦1)-C( 
»2)-C( 
.2)-C( 
.50 
.4 
,8 

ID)+F( 
.1.1)* 
.2)* 
.1)* 

♦  2)* 

♦  D* 
.2)* 
.1)* 
.2)* 


1.1) )*8K.(1.1)  +  (A(2»1 )-C(2»l ) )*BK(2.1) 
1.1) )*BK(1.2)+(A(2»1)-C(2»1 ) )*BK(2»2) 
1.2) )*BK(1,1)+(A(2»2)-C(2.2) )*BK(2»1) 
1.2) )*BK< 1.2)+(A(2.2)-C(2»2) )*B<(2.2) 


IC)*AK( IC-l.ID) 
DE( 1)-E( IB,2.1)*DE(3) 
DE( 1)-E( IB.2,2)*DE(3) 
DE(2)-E( IB,2,1)*DE(4) 
DE(2)-E( IB,2,2)*DE(4) 
DE(5)-E( IB,2,1)*DE(7) 
DE(5)-E( IB,2,2)*DE(7) 
DE(6)-E( IB.2,1)*DE(8) 
DE(6)-E( IB»2,2)*DE(8) 


•  1 

•  1 
.1 
.1 
.1 
.1 
.1 
.8 

STEP( IB)*DP(ID) 

.8 

( ID)+(AK( 1»ID)+2.*AK(2.ID)+2.*AK( 3 ♦ I D)+AK ( 4. ID ) ) /6. 

B) 

(1) 

(2) 

(3) 

(4) 

LV(5) 

LV(6) 

LV(7) 

LV(8) 
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204 


DELK 

TSTE 

TSTE 

THE 

THE 

ELV( 

ELV( 

ELV< 

ELV( 

ELV( 

ELV( 

ELV( 

ELV( 

ELV( 

ELV( 

ELV( 

ELV( 

ELV( 

ELV( 

DO  2 

DO  2 

DO  2 

DE<  I 

DP(1 

DP  (3 

DP  (2 

DP  (4 

ACON 

+  1)> 
ACON 

+2)  ) 
SRT  = 
AKDE 
SDET 
DP  (5 
DP  (6 


IS  THE  MATRIX  OF  DELK  S  AT  T=0 
P( 1)=-TSTEP(1) 
P(2)=-TSTEP(2) 
K  MATRIX*  AS  IS»  IS  AT  T  =  0. 
NEXT  INTEGRATION  IS  FROM  T=0  TO  T=TONE 
1)=BK( 1»1) 
2>»BKU»2) 
3)=BK(2»1) 
4)=BK(2»2) 
5)=0. 
6)=0, 
7)=X0 
8)=Y0 

9)=DELK(1,1) 
10)=DELK(1»2) 
ll)=DELK(2.1) 
12)=DELK(2»2) 
13)=0. 
14)=0. 

02  IA=1,50 

03  IC=1»4 

04  ID=1,14 

D)=ELV( ID)+F( IC)*AK(IC-1,ID) 
)=-A(l,l)*DE( 1)-A(2»1)*DE(3) 
)=-A< 1,2)*DE(1 )-A(2»2)*DE(3) 
>=-A{ 1>1)*DE(2)-A(2.1)*DE(4) 
)=-A( 1,2)*DE(2)-A(2»2)*DE(4) 
(1)=  DE(1)*B(1,1)+DE(3)*B(2»1)+CGUESS*(DE(2)*B( 1 » 1  )+DE ( 4 )*B ( 2 » 


( 2 ) =DE ( 1 ) *B  < 1 1  2 ) +DE ( 3 ) *B ( 2  » 2 ) +  CGUESS*  < DE ( 2 ) *B ( 1 »  2 ) +DE ( 4 ) *B  (  2  » 

SQRTF ( ACON ( 1 ) *ACON ( 1 ) +ACON ( 2 ) *ACON ( 2 ) ) 
T=DE(1)*DE(4)-DE(2)*DE(3) 
=B{ 1,1)*B<2,2)-B(2»1)*B<1,2) 
)=-CGUESS*AKDET*AKDET*BDET*BDET/SRT**3 
)*         AKDET*AKDET*BDET*BDET/SRT**3 


0105 
0106 
0107 
0108 
0109 
0110 
0111 
0112 
0113 
0114 
0115 
0116 
0117 
0118 
0119 
0120 
0121 
0122 
0123 
0124 
0125 
0126 
0127 
0128 
0129 
0130 
0131 
0132 
0133 
0134 
0135 
0136 
0137 
0138 
0139 
0140 
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DP(7)=A(ltl)*DE(7)+A(l»2)*DE(8)+(B(ltl)*ACON( 1 ) +  B ( 1 . 2 ) *ACON ( 2 ) )  /SR    0141 

+T  0142 

DP(8)=A(2»l)*DE(7)+A(2»2)*DE(8)  +  (B(2tD*AC0N(  1  )+B  (  2  »  2  I  *AC0N  (  2  )  )  /SR    0143 

+T  0144 

DP(9)=-A( 1»1)*DE(9)-A(2»1)*DE( 11)  0145 

DP( 11)=-A( 1>2)*DE(9)-A(2»2)*DE( 11)  0146 

DP(  10)=-A(  1,1)*DE(10)-A(2»1)*DF(12)  0147 

DP( 12)=-A( 1,2)*DE( 10)-A(2»2)*DE( 12 )  0148 

DINTFAC=CGUESS*(DE( 1)»DE( 12 ) ~DE ( 4 ) *DE ( 9 ) +DE ( 2 ) *DE ( 1 1 ) -DE ( 3 ) *DE ( 10)    0149 

1      +CGUESS*(DE(2)*DE(12)-DE(4)*DE( 10) ) ) -DE ( 1 ) *DE< 11 )+DE ( 3 ) *DE ( 9 )    0150i 

DP( 13)=-CGUESS*BDET*BDET*AKDET*DINTFAC/SRT**3  0151 

DP(14)=         BDET*BDET*AKDET*DINTFAC/SRT**3  0152 

DO  203  ID=1»14  0153' 

203  AK( IC» ID)=TSTEP(2)*DP( ID)  0154. 

DO  206  ID=1,14  0155  > 

20  6  ELV( ID)=ELV( ID)+(AK( 1, I D)+2.*AK ( 2 » ID )+2 .*AK ( 3 . I D) +AK ( 4 » ID) ) /6.  0156 

BK(ltl)*ELV(l)  0157 

BK( 1.2)=ELV(2)  0158 

BK(2»1)=ELV(3)  0159 

BK(2,2)=ELV<4)  0160 

C0RR(1)=ELV(5)  0161 

CORR(2)=ELV(6)  0162 

XM0D=ELV(7)  0163' 

YM0D=ELV(8)  0l64i 

DELK(1»1)=ELV(9)  0165 

DELK(1»2)=ELV(10)  0166 

DELK(2»1)=ELV( 11)  0167' 

DELK(2.2)=ELV(12)  0168' 

CORR(3)=ELV(13)  0169 

C0RR(4)=ELV(14)  0170 

202  T=T+TSTEP(2)  0171 

CORR  3  AND  4  GIVE  THE  VARIATIONS  OF  X  AND  Y  AT  Tl  DUE  TO  DTI          0172 

GET  K  INVERSE  AT  TONE  AND  F  MINUS  0173 

AKDET=BK(1,1)»BK(2»2)-BK(1,2)»BK(2,1 )  0174. 

XD0T1=DP(7)  0175 

YD0T1=DP(8)  0176 


46 


CKI (1.1 )=BK(2»2)/AKDET  0177 

CKI (1.2)=-BK(2»1)/AKDET  0178 

CKI (2.1)=-BK(1,2)/AKDET  0179 

CKI (2,2)=BK( ltl)/AK0ET  0180 

FLESU  )=ACON(l)/SRT  0181 

FLES(2)=ACON(2)/SRT  0182 

Y0NE(3)=C0RR(4)  0183 

X0NE(3)=CKI (l»l)*CORR(3)+CKI ( 1 ,2 ) *CORR ( 4 )  0184 

XONE(l)=XONE(3)+XDOTl  0185 

XONE(2)=      (CKK l»l)*CORR(l)+CKI ( 1 ,2 ) *CORR ( 2 ) )  0186 

YONE(2)=CORR<2)  0187 

XINMID=XMOD  01*8 

PRINT  559.XINMID  0189 

559  FORMAK20H  559   X  AT  TONE  IS   ,F20.10)  0190 

INTEGRATE  FROM  Tl  TO  CAPT  0191 

ELV(1)=BK(1,1)  0192 

ELV(2)=BK<1»2)  0193 

ELV(3)=BK(2,1)  0194 

ELV(4)=BK(2.2)  0195 

ELV(5)=XMOD  '  0196 

ELV(6)=YMOD  0197 

ELV(7)=0.  0198 

ELV(8)»0.  0199 
DO  302  IA=1,50                                                        •  0200 

DO  303  IC=1,4  0201 

DO  304  ID=1»8  0202 

304  DE( ID)=ELV( ID)+F( IC> *AK ( IC-1 » TO)  0203 

DP(1)=-C(1»1)*DE<1)-C(2#1)*DE(3)  02  04 

DP(3)=-C(1.2)*DE(1 )-C(2»2)*DE(3)  0205 

DP(2)=-C(1»1)*DE(2)-C(2»1)*DE(4)  02  06 

DP(4)=-C(1»2)*DE(2)-C(2»2)*DE{4)  02  07 
AC0N(3)=DE( l)*D(l»l)+DE(3)*D(2»l)+DGUESS*(DE(2)*D(ltl)+DE(4)*D(2tl    02  08 

+))  0209 
AC0N(4)=DE( 1)*D(1»2)+DE(3)*D(2»2)+DGUESS*(DE(2)*D(1.2)+DE(4)*D(2,2    0210 

+))  0211 

SRT=SQRTF(AC0N(3)*AC0N(3)+AC0N(4)*AC0N(4> )  0212 
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CKDET=DE(1)*DE(4)-DE(2)*DE(3)  0213 

DOET=D( ltl)*D(2t2)-D(l»2)*D(2, 1)  0214 
DP(5)=C(l»l)*DE(5)+CQ»2)*DE(6)  +  (D(ltl)*ACON(3)+D(lt2)*ACON(4)  ) /SR    0215 

+  T  0216 
DP(6)=C(2tl>*DE(5>+C(2»2>*DE(6)  +  (D(2»l>*AC0N(3>+D(2»2)*AC0N(4) ) /SR    0217'  ' 

+  T  02181 

DP(7)=         DDET*DDET*CKDET#C<DET/SRT**3  0219' 

DP(8)=-DGUESS*DDET*DDET*CKDET»CKDET/SRT*»3  02  20 

DO  303  ID=1,8  0221 

303  AK( IC» ID)=TSTEP( 1)*DP< ID)  0222" 

IF( IA-1 )801»801»802  0223' 

801  CONTINUE  0224 
XDOT2=DP(5)  0225- 
YD0T2=DP(6)  0226 

802  CONTINUE  0227 
DO  306  ID=1,8  0228 

30  6  ELV( ID)=ELV( ID)+(AK( 1 ♦ ID) +2.*AK ( 2 • I D )+2 .*AK( 3 ♦ I D) +AK ( 4. ID) )/6.  02  29 

302  T=T+TSTEP(1)  0230 

BKA(1,1)=ELV(1)  0231 

BKA(1,2)=ELV(2)  0232 

BKA(2,1)=ELV(3)  0233 

BKA(2,2)=ELV(4)  0234 

XMOD=ELV(5)  0235 

YMOD=ELV(6)  0236 

CORA=ELV(7)  0237 

CORB=ELV(8)  0238 

THIS  COMPLETES  THE  INTEGRATION  0239 

XONE(5)=CORR(3)+BK(ltl)*(XDOTl-XDOT2)+BK(2»l ) ♦< YDOT 1-YDOT2 )  0240 

YONE(l)=CORR(4)+BKUt2)*(XDOTl-XDOT2)+BK(2»2)*(YDOTl-YDOT2)  0241 

XBND=XMOD  0242 

YBND=YMOD  0243 

PRINT  558,XBND,YBND  0244 

558  F0RMAT(4H  558»9H   X(T)«   »F20.10,9H   Y(T)»   .F20.10)  0245 

XDIFF=XMID-XINMID  0246 

YDIFF=YFINAL-YBND  0247 

IF(ABSF(XDIFF)+ABSF(YDIFF)-ERROR)998.998t600l  02  48 


48 


>001  DETER=X0NE< 1)*CORR(2)-XONE(2)*YONE(1)  0249 

IF ( ABSF ( DETER )-l.E-5) 6002 t 6002 .6011  02  50 

Oil  DELTAC=XDIFF/X0NE(2)  0251 

DELTAD=(YDIFF-C0RR(2)*DELTAC)/C0RA  02  52 

DELTAT=0«  0253 

i670  IF(ABSF(DELTAC>-.3)6671»6672»6672  0254 

>672  DELTAC=.5*DELTAC  0255 

DELTAD=.5*DELTAD  0256 

DELTAT=.5*DELTAT  0257 

GO  TO  6670  0258 

>671  GO  TO  9001  0259 

>002  DELTAT  =  XDIFF/XONE(D  0260 

DELTAD=(YDIFF-Y0NE(1)*XDIFF/X0NE(1> ) /CORA  0261 

J673  IF(ABSF(DELTAD)-.3)6674t6675»6675  0262 

>675  DELTAD=.5*DELTAD  0263 

DELTAC=.5*DELTAC  0264 

DELTAT=.5*DELTAT  0265 

GO  TO  6673  0266 

S674  CONTINUE  0267 

DELTAC=0.                                                    .   .  0268 

DETER=XONE( 1)*C0RA  0269 

IF(ABSF(DETER)-1.E-5)9198,9198,9001  0270 

5001  IF(ABSF(DELTAT)-.5)990»9002»9002  0271 

?002  DELTAT=.5*DELTAT               »  0272 

DELTAC=.5*DELTAC  0273 

DELTAD=.5*DEI_TAD  0274 

GO  TO  9001  0275 

990  PRINT  562»DELTAC»DELTAT,DELTAD  0276 

562  FORMAT(12H   DELTAC  =   »F20.10,14H   DELTA  Tl  =   »F20.10.12H   DELTAD    0277 

+  =   •F20.10)  0278 

TTW0=T0NE+DELTAT  0279 

IF(TTWO)996»996»309  02R0 

996  T0NE=T0NE*.5  0281 

CGUESS=CGUESS+(TONE/DELTAT*DELTAC)  02  82 

DGUESS=DGUESS+(TONE/DELTAT*DELTAD)  02  83 

GO  TO  310  0284 
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309  IF(TTWO-CAPT)311,995,995  028!:' 
995  TONE=(TONE+CAPT)*.5  0286 

CGUESS=CGUESS+(CAPT-TONE)*.5/DELTAT*DELTAC   ■  0281 

DGUESS=DGUESS+(CAPT-TONE)*.5/DELTAT*DELTAD  02  8? i: 

GO  TO  310  028* 

311  TONE=TTWO  029C 

CGUESS=CGUESS+DELTAC  029T: 

DGUESS=DGUESS+DELTAD  0292 

310  CONTINUE  029; 
IF (CAPT-TONE-l.E-6)977 1,977 1,9772  02  9*| 

9771  PRINT  9781 

9781  FORMAT(5H  9781tl6H   TONE  TOO  LARGE) 
GO  TO  9991 

9772  IF (TONE-1. £-6)9773,9773,  9774 

9773  PRINT  9782 

9782  FORMAT* 5H  9782, 16H   TONE  TOO  SMALL) 
GO  TO  9991 

9774  CONTINUE 
322  CONTINUE 

PRINT  563, CGUESS, TONE 
563  FORMAT  ( 16H  563   NEW  C  IS   ,F20.10,15H   NEW  TONE  IS   ,F20.10) 
999  CONTINUE 

GO  TO  1 
99  8  DTMT=XONE(2)*(CORA*XONE(5)-CORB*YONE( 1 ) ) 
+     +XONE( 1 )*(CORB*CORR(2)-CORA*CORR(l ) ) 

IF (ABSF(DTMT)-. 01) 960 1,960 1,9602 

9601  IF(ABSF(DTMT)-. 00001)9603, 9603, 9604 

9603  PRINT  9651, XBND 
9651  FORMAT(5H  9651, 9H   XMAX*   ,F20.10) 

GO  TO  9991 

9604  CONTINUE 
GO  TO  9602 

9602  EPS2=EPSFAC*DTMT 
DELTAX=XBND-X11 
IF ( DELTAX ) 8401 , 8401 , 8402 

8401  EPS=EPS*.l 
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GO  TO  9627  0321 

^402  EPS=EPS2  0322 

X11=XBND  0323 

GO  TO  9616  0324 

?627  T0NE=T0NE-DELT1T  0325 

DGUESS=DGUESS-DELT1D  0326 

CGUESS=CGUESS-DELT1C  0327 

9616  DELTlOEPS*(-X0NE(l>*C0RA)  0328 
DELT10=EPS*(XONE( 1 ) *CORR( 2 )-X0NE ( 2 ) *Y0NE < 1 ) )  03  29 
DELT1T=EPS*(X0NE(2 )*CORA)  0330 

621  IF (ABSF(DELT1C)+ABSF(DELT1D)+ABSF(DELT1T)-. 5) 6622 .6623*662  3  03  31 

6623  DELT1C=.5*DELT1C  0332 

DELT1D=.5*DELT1D  0333 

DELT1T=.5*DELT1T  0334 

GO  TO  6621  0335 

6622  CONTINUE  0336 

PRINT  9611.DELT1C.DELT1D.DELT1T  0337 

9611  FORMAT(5H  9601»10H  DELTAC"   .F20.10.11H   DELTAD=   .F20.10.11H   DEL    0338 

1TAT=   .F20.10)  0339 

CGUESS=CGUESS+DELT1C  0340 

DGUESS=DGUESS+DELT1D  0341 

T0NE=T0NE+DELT1T  0342 

ATANC=ATANF(CGUESS)  0343 

ATAND=ATANF(DGUESS)  0344 

PRINT  9617»ATANC.ATAND  0345 

9617  FORMAT(5H  9617»6H   C=   »F20*lO»10H  RADIANS.  »6H   D=   .F20.10.9H   R    0346 
|    +ADIANS)  0347 

DELTC=ATANC-ATC11  0348 

DELTD=ATAND-ATD11  0349 

IF(ABSF(DELTC)+ABSF(DELTD)-l.E-6)  9981,9981.9982                     0350 

9981  PRINT  9983  0351 

9983  FORMAT ( 5H  9983.49H   CORRECTION  IS  LESS  THAN  .000001    -  WE  ARE  DON    0352 

IE//)  0353 

GO  TO  9991  0354 

9982-ATCll=ATANC  0355 

ATD11=ATAND  0356 
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9985    CONTINUE  03^7 

PRINT    9612, EPS  03^8 

9612    FORMAT(12H    9612    EPS*       ,F20.10)  0359 

AAK ( 1 , 1 ) =BK ( 1 ♦ 1 )+CGUESS*BK (1,2)  0360 

AAK(2,1)=BK(2,1)+CGUESS*BK(2,2)  0361 

AAK ( 1 ,2 ) =BK ( 1 .1 )+DGUESS*BK (1,2)  0362 

AAK  <  2 ,2 ) =BK ( 2 .1 )+DGUESS*BK ( 2  t2 )  0363 

PRINT    9987  0364 

9987  FORMAT(5H  9987, 18H  K*  AT  TONE  MINUS* 10X t 15HK*  AT  TONE  PLUS//)  0365 
PRINT    9988, ( (AAK( I .J) , J=l ,2 ) » I *1 ,2 )  0366 

9988  FORMAT( 10X, F12. 5, 14X ,F12. 5 )  0367 
PRINT    9989  0368 

9989  FORMAT(1HO,12H  XOOT  MINUS, 10x ,9HXDOT  PLUS)  0369 
PRINT    9990,XDOT1,XDOT2,YDOT1,YDOT2  0370 

9990  FORMAT(5X,F12.5,10X,F12»5)  0371 
GO    TO    1  0372 

9198  PRINT    9199  0373 

9199  FORMAT(68H    RANK   OF    FIRST    CORRECTION    MATRIX    IS    ONE    -    NO    ADMISSIBLE  0374 
-      1CURVES    EXIST)  0375 

9991  CONTINUE  0376 
STOP  0377 
END  0378 
END  0379 
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